#include "device.h"

/* --------------------------------------------- */
/* --------------------------------------------- */

void Device::set_strain(short flag) {
  
  int k,i,j;
  double d,a_temp,tmp;

  double RELAX_LENGTH = sim->relax_length;

  if(lyr[0].agrid->a(0,0)==0) return;

  if(flag) {
    if(ynum_mesh_points<=1) {
      for(i=0;i<lyr[0].numx();i++) {
	lyr[0].grid1d[i].set_a(lyr[0].agrid->a(lyr[0].grid1d[i].ret_alloy_x(),lyr[0].grid1d[i].ret_alloy_y()));
	lyr[0].grid1d[i].set_c(lyr[0].agrid->c(lyr[0].grid1d[i].ret_alloy_x(),lyr[0].grid1d[i].ret_alloy_y()));
      }
      for(k=1;k<num_lyr;k++) {
	a_temp = lyr[k-1].grid1d[lyr[k-1].numx()-1].ret_a();
	if( lyr[k].agrid->a(lyr[k].grid1d[0].ret_alloy_x(),lyr[k].grid1d[0].ret_alloy_y())!= a_temp ) {
	  d = 0;
	  for(i=0;i<lyr[k].numx();i++) {
	    tmp = lyr[k].agrid->a(lyr[k].grid1d[i].ret_alloy_x(),lyr[k].grid1d[i].ret_alloy_y());
	    if(d<RELAX_LENGTH) lyr[k].grid1d[i].set_a(a_temp);
	    else lyr[k].grid1d[i].set_a(tmp + (a_temp-tmp)*RLX(d,fabs(a_temp-tmp)/tmp));
	    lyr[k].grid1d[i].set_c(lyr[k].grid1d[i].ret_a(),lyr[k].grid1d[i].ret_alloy_x(),lyr[k].grid1d[i].ret_alloy_y());
	    d += lyr[k].grid1d[i].xsize();
	  }
	}
	else for(i=0;i<lyr[k].numx();i++) {
	  lyr[k].grid1d[i].set_a(lyr[k].agrid->a(lyr[k].grid1d[i].ret_alloy_x(),lyr[k].grid1d[i].ret_alloy_y()));
	  lyr[k].grid1d[i].set_c(lyr[k].agrid->c(lyr[k].grid1d[i].ret_alloy_x(),lyr[k].grid1d[i].ret_alloy_y()));
	}
      }
    }
    else if(ynum_mesh_points>1) {
      for(i=0;i<lyr[0].numx();i++) {
	for(j=0;j<lyr[0].numy();j++) {
	  lyr[0].grid2d[i][j].set_a(lyr[0].agrid->a(lyr[0].grid2d[i][j].ret_alloy_x(),lyr[0].grid2d[i][j].ret_alloy_y()));
	  lyr[0].grid2d[i][j].set_c(lyr[0].agrid->c(lyr[0].grid2d[i][j].ret_alloy_x(),lyr[0].grid2d[i][j].ret_alloy_y()));
	}
      }
      for(k=1;k<num_lyr;k++) {
	for(j=0;j<lyr[0].numy();j++) {
	  a_temp = lyr[k-1].grid2d[lyr[k-1].numx()-1][j].ret_a();
	  if( lyr[k].agrid->a(lyr[k].grid2d[0][j].ret_alloy_x(),lyr[k].grid2d[0][j].ret_alloy_y())!= a_temp ) {
	    d = 0;
	    for(i=0;i<lyr[k].numx();i++) {
	      tmp = lyr[k].agrid->a(lyr[k].grid2d[i][j].ret_alloy_x(),lyr[k].grid2d[i][j].ret_alloy_y());
	      if(d<RELAX_LENGTH) lyr[k].grid2d[i][j].set_a(a_temp);
	      else lyr[k].grid2d[i][j].set_a(tmp + (a_temp-tmp)*RLX(d,fabs(a_temp-tmp)/tmp));
	      lyr[k].grid2d[i][j].set_c(lyr[k].grid2d[i][j].ret_a(),lyr[k].grid2d[i][j].ret_alloy_x(),
					lyr[k].grid2d[i][j].ret_alloy_y());
	      d += lyr[k].grid2d[i][j].xsize();
	    }
	  }
	  else { 
	    for(i=0;i<lyr[k].numx();i++) {
	      lyr[k].grid2d[i][j].set_a(lyr[k].agrid->a(lyr[k].grid2d[i][j].ret_alloy_x(),lyr[k].grid2d[i][j].ret_alloy_y()));
	      lyr[k].grid2d[i][j].set_c(lyr[k].agrid->c(lyr[k].grid2d[i][j].ret_alloy_x(),lyr[k].grid2d[i][j].ret_alloy_y()));
	    }
	  }
	}
      }
    }
  }
  else {
    if(ynum_mesh_points<=1) {
      for(k=0;k<num_lyr;k++) 
	for(i=0;i<lyr[k].numx();i++) {
	  lyr[k].grid1d[i].set_a(lyr[k].agrid->a(lyr[k].grid1d[i].ret_alloy_x(),lyr[k].grid1d[i].ret_alloy_y()));
	  lyr[k].grid1d[i].set_c(lyr[k].agrid->c(lyr[k].grid1d[i].ret_alloy_x(),lyr[k].grid1d[i].ret_alloy_y()));
	}
    }
    else if(ynum_mesh_points>1) {
      for(k=0;k<num_lyr;k++) 
	for(i=0;i<lyr[k].numx();i++) 
	  for(j=0;j<lyr[k].numy();j++) {
	    lyr[k].grid2d[i][j].set_a(lyr[k].agrid->a(lyr[k].grid2d[i][j].ret_alloy_x(),lyr[k].grid2d[i][j].ret_alloy_y()));
	    lyr[k].grid2d[i][j].set_c(lyr[k].agrid->c(lyr[k].grid2d[i][j].ret_alloy_x(),lyr[k].grid2d[i][j].ret_alloy_y()));
	  }
    }
  }
}

/* ---------------------------------------------------------------- */
/* ---------------------------------------------------------------- */
Device::Device(double PZ, double PSP) {

  int * numx, * numy, * numt;
  double * length,* width,* time, * doping_lhs,* doping_rhs,* perx_lhs,* perx_rhs,* pery_lhs,* pery_rhs,* bias;

  FILE *fp;
  int i;

  /* INITIALISE CLASS CONDITION */
  sim = new Condition ("INPUT_FILES/simulation.inp");
  if(PZ!=-1) sim->Pz_percentage = PZ;
  if(PSP!=-1) sim->Psp_percentage = PSP;
  sim->show();
  if(sim==NULL) nrerror("Device: memory allocation error for sim");

  /* READ THE DEVICE.INP FILE */
  if ( (fopen_s(&fp,"INPUT_FILES/device.inp","r"))==NULL ) nrerror("READER_DEV: Cannot open file");
  skip_line(fp);
  fscanf_s(fp,"%d",&num_lyr);

  /* Allocate temporary arrays */
  length = vector(1,num_lyr);
  width = vector(1,num_lyr);
  time = vector(1,num_lyr);
  doping_lhs = vector(1,num_lyr);
  doping_rhs = vector(1,num_lyr);
  perx_lhs = vector(1,num_lyr);
  perx_rhs = vector(1,num_lyr);
  pery_lhs = vector(1,num_lyr);
  pery_rhs = vector(1,num_lyr);
  numx = ivector(1,num_lyr);
  numy = ivector(1,num_lyr);
  numt = ivector(1,num_lyr);
  
  skip_line(fp);
  skip_line(fp);
  printf("\n\n");
  for(i=1;i<=num_lyr;i++) { 
    fscanf_s(fp,"%d%d%d%le%le%le%le%le%lf%lf%lf%lf"
	   ,&numx[i],&numy[i],&numt[i],&length[i],&width[i],&time[i],&doping_lhs[i],&doping_rhs[i]
	   ,&perx_lhs[i],&perx_rhs[i],&pery_lhs[i],&pery_rhs[i]);
    printf("%d->(%d,%d,%d)  x-lgth=%7.2le  y-lgth=%7.2le  t-lgth=%7.2le  dop=(%7.2le,%7.2le)  %%=(%5.2lf,%5.2lf)  %%=(%5.2lf,%5.2lf)\n"
	   ,i,numx[i],numy[i],numt[i],length[i],width[i],time[i],doping_lhs[i],doping_rhs[i]
	   ,perx_lhs[i],perx_rhs[i],pery_lhs[i],pery_rhs[i]);
  }
  printf("\n\n");
  fclose(fp);

  /* DISTRIBUTE THE BIAS ON THE LAYERS  */
  bias = vector(1,num_lyr);
  bias[1] = - 0.5 * sim->initial_bias;
  bias[num_lyr] =  0.5 * sim->initial_bias;
  for(i=2;i<num_lyr;i++) {
    if(doping_rhs[1]>0 && doping_lhs[i]>0) {
      if(doping_lhs[i]*1e2<doping_rhs[1]) bias[i] = 0;
      else bias[i] = bias[1];
    }
    else if(doping_rhs[1]<0 && doping_lhs[i]<0) {
      if(-doping_lhs[i]*1e2<-doping_rhs[1]) bias[i] = 0;
      else bias[i] = bias[1];
    }
    else if(doping_lhs[num_lyr]<0 && doping_rhs[i]<0) {
      if(-doping_rhs[i]*1e2<-doping_lhs[num_lyr]) bias[i] = 0;
      else bias[i] = bias[num_lyr];
    }
    else if(doping_lhs[num_lyr]>0 && doping_rhs[i]>0) {
      if(doping_rhs[i]*1e2<doping_lhs[num_lyr]) bias[i] = 0;
      else bias[i] = bias[num_lyr];
    }
    else {
      nrerror("Something wrong with the bias setting!");
    }
  }

  /* CREATE AND INITIALISE THE LAYERS  */
  lyr = new Layer [num_lyr];
  if(lyr==NULL) nrerror("Device: memory allocation error for lyr. OUT OF MEMORY");
  lyr[0].set(numx[1],numy[1],numt[1],length[1],width[1],time[1],doping_lhs[1],doping_rhs[1]
	     ,perx_lhs[1],perx_rhs[1],pery_lhs[1],pery_rhs[1],bias[1]
	     ,sim->traps,sim->barrier_left_contact,sim->input_file_name,sim->Pz_percentage,sim->Psp_percentage
	     ,sim->fluctuation,sim->fluctuation_percentage);
  for(i=2;i<num_lyr;i++) 
    lyr[i-1].set(numx[i],numy[i],numt[i],length[i],width[i],time[i],doping_lhs[i],doping_rhs[i],
		 perx_lhs[i],perx_rhs[i],pery_lhs[i],pery_rhs[i],bias[i],sim->traps,0,sim->input_file_name,sim->Pz_percentage
		 ,sim->Psp_percentage,sim->fluctuation,sim->fluctuation_percentage);
  lyr[num_lyr-1].set(numx[num_lyr],numy[num_lyr],numt[num_lyr],length[num_lyr],width[num_lyr],time[num_lyr],doping_lhs[num_lyr],
		     doping_rhs[num_lyr],perx_lhs[num_lyr],perx_rhs[num_lyr],pery_lhs[num_lyr],pery_rhs[num_lyr],
		     bias[num_lyr],sim->traps,sim->barrier_right_contact,
		     sim->input_file_name,sim->Pz_percentage,sim->Psp_percentage,sim->fluctuation,sim->fluctuation_percentage);

  /* COUNT THE TOTAL NUMBER OF MESH POINTS  */
  xnum_mesh_points = 0;
  for(i=1;i<=num_lyr;i++) xnum_mesh_points += numx[i];
  ynum_mesh_points = numy[1];
  tnum_mesh_points = numt[1];
  
  /* SET THE STRAIN */
  set_strain(1);
  printf("Strain set\n");

  /* FREE ALL THE TEMPORARY ARRAYS */
  free_vector(length,1,num_lyr);
  free_vector(width,1,num_lyr);
  free_vector(time,1,num_lyr);
  free_vector(doping_lhs,1,num_lyr);
  free_vector(doping_rhs,1,num_lyr);
  free_vector(perx_lhs,1,num_lyr);
  free_vector(perx_rhs,1,num_lyr);
  free_vector(pery_lhs,1,num_lyr);
  free_vector(pery_rhs,1,num_lyr);
  free_vector(bias,1,num_lyr);
  free_ivector(numx,1,num_lyr);
  free_ivector(numy,1,num_lyr);
  free_ivector(numt,1,num_lyr);


}


/* ---------------------------------------------------------------- */
/* Write results to output files
/* ---------------------------------------------------------------- */

void Device::dump_files(char ch) {

  int i;
  double pos;
  char namefile[9][25]; // array of file names
  FILE * fp1, *fp2, *fp3, *fp4, *fp5 , *fp6, *fp7, *fp8, *fp9;

  strcpy_s(namefile[0],"OUTPUT_FILES/bandx.dat");
  strcpy_s(namefile[1],"OUTPUT_FILES/carrx.dat");
  strcpy_s(namefile[2],"OUTPUT_FILES/potex.dat");
  strcpy_s(namefile[3],"OUTPUT_FILES/currx.dat");
  strcpy_s(namefile[4],"OUTPUT_FILES/qwlex.dat");
  strcpy_s(namefile[5],"OUTPUT_FILES/qwlhx.dat");
  strcpy_s(namefile[6],"OUTPUT_FILES/levex.dat");
  strcpy_s(namefile[7],"OUTPUT_FILES/levhx.dat");
  strcpy_s(namefile[8],"OUTPUT_FILES/optix.dat");

  for(i=0;i<9;i++) namefile[i][17]=ch;

  fopen_s(&fp1,namefile[0],"w");
  fopen_s(&fp2,namefile[1],"w");
  fopen_s(&fp3,namefile[2],"w");
  fopen_s(&fp4,namefile[3],"w");
  fopen_s(&fp5,namefile[4],"w");
  fopen_s(&fp6,namefile[5],"w");
  fopen_s(&fp7,namefile[6],"w");
  fopen_s(&fp8,namefile[7],"w");
  fopen_s(&fp9,namefile[8],"w");

  if( !fp1 || !fp2 || !fp3 || !fp4 || !fp5 || !fp6 || !fp7 || !fp8|| !fp9) {
    printf("DEVICE::DUMP_FILES I cannot open one of the files\n");
    exit(1);
  }
  printf("\nDEVICE::DUMP_FILES ...\n");
  pos = 0;
  for(i=0;i<num_lyr;i++) {
    lyr[i].show(pos,fp1,fp2,fp3,fp4,fp5,fp6,fp7,fp8,fp9);
    pos += lyr[i].xsize();
  }
  fclose(fp1);
  fclose(fp2);
  fclose(fp3);
  fclose(fp4);
  fclose(fp5);
  fclose(fp6);
  fclose(fp7);
  fclose(fp8);
  fclose(fp9);

	
}


/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */

void Device::Poisson_solvde(short debug) {
  
  int 	i,*indexv,NUM;
  const int ITMAX	=	5000;		/* maximum number of iterations */
  const double conv	=	1.0e-11;	/* convergence criterion	*/
  double slowc = 0.01;
  double * scalv,**s,***c,**y;

  NUM = xnum_mesh_points-2; /* Skip the first and the last mesh point */

  /* ALLOCATE AND SET DYNAMIC ARRAYS */

  if(debug) printf("POISSON SOLVER ...\n");

  scalv = vector(1,NE);	/*Contains average values for dependent variables. It is used to scale the error */
  scalv[1] = 1.5/VT;
  scalv[2] = (1.0e-4)*L*L/Q;

  indexv = ivector(1,NE); /*Can be used to switch the order of the solved variables */
  indexv[1] = 1;
  indexv[2] = 2;

  y = matrix(1,NE,1,NUM);    /* Dependent variable vector */
  s = matrix(1,NE,1,2*NUM+1);
  c = tensor(1,NE,1,NE-NB+1,1,NUM+1);

  /*	SET INITIAL CONDITIONS */

  for(i=1;i<=NUM;i++) {
    y[1][i] = phi(i)/VT;
    y[2][i] = D(i)*L*L/Q;
  }

  /* SOLVDE */

  solvde(ITMAX,conv,slowc,scalv,indexv,NE,NB,NUM,y,c,s,*this);

  /* UPDATE POTENTIAL AND DISPLACEMENT IN THE MESH */
  /* --------------------------------------------- */

  for(i=1;i<=NUM;i++) {
    set_phi(i,y[1][i]*VT);
    set_D(i,y[2][i]*Q/(L*L));
  }
 

  /* DEALLOCATE MEMORY */
  /* ----------------- */
  free_vector(scalv,1,NE);
  free_ivector(indexv,1,NE);
  free_matrix(y,1,NE,1,NUM);
  free_matrix(s,1,NE,1,2*NUM+1);
  free_tensor(c,1,NE,1,NE-NB+1,1,NUM+1);

  if(debug) printf("... DONE!\n\n");

}

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */


void Device::difeq( int k,int k1,int k2,int jsf,int is1,int isf,int indexv[],int ne,double **s, double **y ) {

  int i,j;
  double aa0,aa1,b0,b1,c0,c1,xx;
  Phys_var * p0, *p1;


  if(k==k1) {
    s[2][ne+indexv[1]] = 1.0;
    s[2][ne+indexv[2]] = 0.0;
    s[2][jsf] =	y[1][1]-(phi(1)/VT);
  }
  else if(k>k2) {
    s[1][ne+indexv[1]] = 1.0;
    s[1][ne+indexv[2]] = 0.0;
    s[1][jsf] = y[1][k2]-(phi(xnum_mesh_points)/VT);
  }
  else {
    MeshToLayer(k-1,i,j);
    p0 = & lyr[i].grid1d[j];
    MeshToLayer(k,i,j);
    p1 = & lyr[i].grid1d[j];
    //aa1 = p1->nqw()+p1->pqw();
    //aa0 = p0->nqw()+p0->pqw();
    aa1 = 0;
    aa0 = 0;
    b1 = -((Q/EPS0)/(L*VT))*(1.0/p1->eps());
    b0 = -((Q/EPS0)/(L*VT))*(1.0/p0->eps());
    c1 = (L/(VT*EPS0))*(p1->Polarization()/p1->eps());
    c0 = (L/(VT*EPS0))*(p0->Polarization()/p0->eps());
    xx = 0.5*(p1->xsize()+p0->xsize())/L;

    if(sim->traps) {
      s[2][indexv[1]] = -VT*xx*0.5*L*L*L*(p0->dND(VT*y[1][k-1])-p0->dNA(VT*y[1][k-1])
					  -p0->dntrap(VT*y[1][k-1])-p0->dn(VT*y[1][k-1])+p0->dp(VT*y[1][k-1]) );
      s[2][indexv[2]] = -1.0;
      s[2][ne+indexv[1]] = -VT*xx*0.5*L*L*L*(p1->dND(VT*y[1][k])-p1->dNA(VT*y[1][k])
					     -p1->dntrap(VT*y[1][k])-p1->dn(VT*y[1][k])+p1->dp(VT*y[1][k]) );
      s[2][ne+indexv[2]] = 1.0;      
      s[2][jsf] = (y[2][k]-y[2][k-1]) -xx*0.5*L*L*L*( aa1+aa0+p0->ND(VT*y[1][k-1])+p1->ND(VT*y[1][k])-p0->NA(VT*y[1][k-1])-p1->NA(VT*y[1][k])
						      -p0->n(VT*y[1][k-1])-p1->n(VT*y[1][k])+p0->p(VT*y[1][k-1])+p1->p(VT*y[1][k])
						      -p0->ntrap(VT*y[1][k-1])-p1->ntrap(VT*y[1][k]));
    }
    else {
      s[2][indexv[1]] = -VT*xx*0.5*L*L*L*(p0->dND(VT*y[1][k-1])-p0->dNA(VT*y[1][k-1])-p0->dn(VT*y[1][k-1])+p0->dp(VT*y[1][k-1]) );
      s[2][indexv[2]] = -1.0;
      s[2][ne+indexv[1]] = -VT*xx*0.5*L*L*L*(p1->dND(VT*y[1][k])-p1->dNA(VT*y[1][k])-p1->dn(VT*y[1][k])+p1->dp(VT*y[1][k]) );
      s[2][ne+indexv[2]] = 1.0;
      s[2][jsf] = (y[2][k]-y[2][k-1]) -xx*0.5*L*L*L*( aa1+aa0+p0->ND(VT*y[1][k-1])+p1->ND(VT*y[1][k])-p0->NA(VT*y[1][k-1])-p1->NA(VT*y[1][k])
						      -p0->n(VT*y[1][k-1])-p1->n(VT*y[1][k])+p0->p(VT*y[1][k-1])+p1->p(VT*y[1][k]));
    }

    s[1][indexv[1]] = -1.0;
    s[1][indexv[2]] = -xx*0.5*b0;
    s[1][ne+indexv[1]] = 1.0;
    s[1][ne+indexv[2]] = -xx*0.5*b1;
    s[1][jsf] =	(y[1][k]-y[1][k-1])-xx*0.5*( b1*y[2][k] + b0*y[2][k-1] + c1 + c0 );
  }
	
}

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */

void Device::mesh_modifier() {
  int k,i,j;
  double d,pot1,pot2,lgth,delta;
  double max_length = 1.0e-5;
  double min_length = 1.0e-8;

  k = 0;
  i = 1;
  lgth = 0;
  while(k<num_lyr) {
    pot1 = lyr[k].grid1d[i-1].ret_phi();
    pot2 = lyr[k].grid1d[i].ret_phi();
    delta = fabs(pot1-pot2);
    if(delta<1e-4) {
      d = lyr[k].grid1d[i].xsize();
      if(d<max_length) {
	lyr[k].grid1d[i].double_xsize();
	lgth += d;
      }
    }
    else if(delta>1e-2) {
      d = lyr[k].grid1d[i].xsize();
      if(d>min_length) {
	lyr[k].grid1d[i].half_xsize();
	lgth += -0.5*d;
      }
    }
    if(i<lyr[k].numx()-1) i++;
    else {
      if(lgth!=0)
	for(j=0;j<lyr[k].numx();j++) lyr[k].grid1d[j].incr_xsize(-lgth/double(lyr[k].numx()));
      k++; 
      i=0;
      lgth = 0;
    }
  }
  
}

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */

void Device::Poisson_2D() {
  
  int i,j,k,iter;
  double PRECISION = 1e-5;
  double lambda ,epsx,epsy,old;
  double dPx, dPy, z, dz;
  double ** A, ** B, ** C, ** D, ** E, ** pot;
  
  printf("Poisson\n");
  pot = matrix(1,xnum_mesh_points,1,ynum_mesh_points);
  A = matrix(2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  B = matrix(2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  C = matrix(2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  D = matrix(2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  E = matrix(2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  
  lambda = Q/(EPS0);

  /* Temporary Copy for working on */
  for(j=1;j<=ynum_mesh_points;j++)
    for(i=1;i<=xnum_mesh_points;i++) {
      pot[i][j] = phi2D(i,j);
    }
      
  for(j=2;j<ynum_mesh_points;j++) {
    for(i=2;i<xnum_mesh_points;i++) {
      epsx = (eps(i+1,j)-eps(i,j))/(0.5*(dx(i+1,j)+dx(i,j)))+((eps(i,j)-eps(i-1,j))/(0.5*(dx(i-1,j)+dx(i,j))));
      epsy = (eps(i,j+1)-eps(i,j))/(0.5*(dx(i,j+1)+dx(i,j)))+((eps(i,j)-eps(i,j-1))/(0.5*(dx(i,j-1)+dx(i,j))));
      A[i][j] = (0.25*epsy/(-dy(i,j-1)-dy(i,j)))+2.0/(0.5*(dy(i,j-1)+dy(i,j))*(0.5*dy(i,j-1)+0.5*dy(i,j+1)+dy(i,j)));
      B[i][j] = (0.25*epsx/(-dx(i-1,j)-dx(i,j)))+2.0/(0.5*(dx(i-1,j)+dx(i,j))*(0.5*dx(i-1,j)+0.5*dx(i+1,j)+dx(i,j)));
      C[i][j] = (0.25*epsx/(-dx(i+1,j)-dx(i,j)))+2.0/(0.5*(dx(i+1,j)+dx(i,j))*(0.5*dx(i-1,j)+0.5*dx(i+1,j)+dx(i,j)));
      D[i][j] = (0.25*epsy/(-dy(i,j+1)-dy(i,j)))+2.0/(0.5*(dy(i,j+1)+dy(i,j))*(0.5*dy(i,j-1)+0.5*dy(i,j+1)+dy(i,j)));
      dPx = ((Polar(i,j)-Polar(i-1,j))/(dx(i,j)+dx(i-1,j)))+((Polar(i+1,j)-Polar(i,j))/(dx(i+1,j)+dx(i,j)));
      dPy = ((Polar(i,j)-Polar(i,j-1))/(dy(i,j)+dy(i,j-1)))+((Polar(i,j+1)-Polar(i,j))/(dy(i,j+1)+dy(i,j)));
      E[i][j] = (dPx + dPy)/EPS0;
    }
  }
   
  for(k=0;k<3;k++) {
    /* Newton Method*/
    for(j=2;j<ynum_mesh_points;j++)
      for(i=2;i<xnum_mesh_points;i++) {
	iter = 0;
	do {
	  z = A[i][j]*(pot[i][j-1]-pot[i][j])
	    + B[i][j]*(pot[i-1][j]-pot[i][j])
	    + C[i][j]*(pot[i+1][j]-pot[i][j])
	    + D[i][j]*(pot[i][j+1]-pot[i][j])
	    - E[i][j] + (Q/EPS0) * charge(i,j,pot[i][j]);
	  dz = - A[i][j]-B[i][j]-C[i][j]-D[i][j]+(Q/EPS0)*dcharge(i,j,pot[i][j]);
	  old = pot[i][j];
	  pot[i][j]=-z/dz;
	  iter++;
	  // printf("Iter = %d  error = %le \n",iter,fabs(pot[i][j]-old));
	} while(fabs(pot[i][j]-old)>PRECISION && iter<1000);
	if(iter>=1000) nrerror("2DPoisson: Precision not reached");
      }
    for(i=1;i<=xnum_mesh_points;i++) {
      pot[i][1]=pot[i][3];
      pot[i][ynum_mesh_points]=pot[i][ynum_mesh_points-2];
    }
  }
  
  /* Update potential */
  for(j=1;j<=ynum_mesh_points;j++)
    for(i=1;i<=xnum_mesh_points;i++) set_phi(i,j,pot[i][j]);
  
  /* Free Matrices */
  free_matrix(A,2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  free_matrix(B,2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  free_matrix(C,2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  free_matrix(D,2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  free_matrix(E,2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  free_matrix(pot,1,xnum_mesh_points,1,ynum_mesh_points);

}


// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

// This function evaluates the net charge that is sum of the 
// ionised dopants, of the free carriers and of the confined and trapped carriers.
// The required parameters are the mesh position x,y and the 
// value of the potential at that mesh point.
// The confined carriers are not considered varying with the potential.

double Device::charge(int numx, int numy, double pot) {

    int k,i,j;
    double temp=0;

    MeshToLayer(numx,k,i);
    j = numy-1;

    temp = lyr[k].grid2d[i][j].ND(pot) - lyr[k].grid2d[i][j].NA(pot)
      + lyr[k].grid2d[i][j].p(pot) - lyr[k].grid2d[i][j].n(pot)
      + lyr[k].grid2d[i][j].ptrap(pot) - lyr[k].grid2d[i][j].ntrap(pot)
	  + lyr[k].grid2d[i][j].pqw() - lyr[k].grid2d[i][j].nqw();

    return temp;

  }

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */

// This is the derivative of the function charge(int,int,double) 
// in the potential. See the function charge(int,int,double).

double Device::dcharge(int numx, int numy, double pot) {

    int k,i,j;
    double temp=0;

    MeshToLayer(numx,k,i);
    j  = numy - 1;

    temp = lyr[k].grid2d[i][j].dND(pot) - lyr[k].grid2d[i][j].dNA(pot)
         + lyr[k].grid2d[i][j].dp(pot) -  lyr[k].grid2d[i][j].dn(pot)
         + lyr[k].grid2d[i][j].dptrap(pot) -  lyr[k].grid2d[i][j].dntrap(pot);

    return temp;

  }

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */

void Device::Continuity_2D() {

  int i,j,jj,k,l,num_x,num_y,iter;
  double **A, **B, **C, **D, **En, **ef;
  double **Ep, **hf, old,tmp1,tmp2,z,dz,x,PRECISION;
  double vr = sqrt((Q/MASS)*VT/PI), gamma = 1.0;
  FILE * fp;

  PRECISION = 1e-5;

  printf("Continuity_2D \n");

  A = matrix(2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  B = matrix(2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  C = matrix(2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  D = matrix(2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  En = matrix(1,xnum_mesh_points,1,ynum_mesh_points);
  Ep = matrix(1,xnum_mesh_points,1,ynum_mesh_points);

  ef = matrix(1,xnum_mesh_points,1,ynum_mesh_points);
  hf = matrix(1,xnum_mesh_points,1,ynum_mesh_points);

  for(j=1;j<=ynum_mesh_points;j++)
    for(i=1;i<=xnum_mesh_points;i++) {
      ef[i][j] = phin(i,j);
      hf[i][j] = phip(i,j);
    }

  for(j=2;j<ynum_mesh_points;j++) {
    for(i=2;i<xnum_mesh_points;i++) {
      A[i][j] = 2.0/(0.5*(dy(i,j-1)+dy(i,j))*(0.5*dy(i,j-1)+0.5*dy(i,j+1)+dy(i,j)));
      B[i][j] = 2.0/(0.5*(dx(i-1,j)+dx(i,j))*(0.5*dx(i-1,j)+0.5*dx(i+1,j)+dx(i,j)));
      C[i][j] = 2.0/(0.5*(dx(i+1,j)+dx(i,j))*(0.5*dx(i-1,j)+0.5*dx(i+1,j)+dx(i,j)));
      D[i][j] = 2.0/(0.5*(dy(i,j+1)+dy(i,j))*(0.5*dy(i,j-1)+0.5*dy(i,j+1)+dy(i,j)));
      En[i][j] = mobn(i,j);
      Ep[i][j] = mobp(i,j);
    }
  }


  for(j=2;j<ynum_mesh_points;j++) {
    jj = j-1;
    for(l=2;l<xnum_mesh_points;l++) {
      iter = 0;
      MeshToLayer(l,k,i);
      printf("l=%d  j=%d \n",l,j);
      do {
	if((k>0 && k<num_lyr-1 && i>0 && i<lyr[k].numx()-1) || (k==0 && i>0 && i<lyr[k].numx()-1 ) || (k==num_lyr-1 && i>0 && i<lyr[k].numx()-1 ) ) {
	  printf("Far from the interfaces\n");
	  z = A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*(nf(l,j-1,ef[l][j-1])+nf(l,j,ef[l][j]))*(-ef[l][j-1]+ef[l][j])
	    + B[l][j]*0.5*(En[l-1][j]+En[l][j])*0.5*(nf(l-1,j,ef[l-1][j])+nf(l,j,ef[l][j]))*(-ef[l-1][j]+ef[l][j])
	    + C[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*(nf(l+1,j,ef[l+1][j])+nf(l,j,ef[l][j]))*(-ef[l+1][j]+ef[l][j])
	    + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*(nf(l,j+1,ef[l][j+1])+nf(l,j,ef[l][j]))*(-ef[l][j+1]+ef[l][j])
	    + lyr[k].grid2d[i][j-1].SRH_Rec_ef(ef[l][j]);
	  dz = A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*dnf(l,j,ef[l][j])*(-ef[l][j-1]+ef[l][j])
	    + A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*(nf(l,j-1,ef[l][j-1])+nf(l,j,ef[l][j]))*ef[l][j]
	    + B[l][j]*0.5*(En[l-1][j]+En[l][j])*0.5*dnf(l,j,ef[l][j])*(-ef[l-1][j]+ef[l][j])
	    + B[l][j]*0.5*(En[l-1][j]+En[l][j])*0.5*(nf(l-1,j,ef[l-1][j])+nf(l,j,ef[l][j]))*ef[l][j]
	    + C[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*dnf(l,j,ef[l][j])*(-ef[l+1][j]+ef[l][j])
	    + C[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*(nf(l+1,j,ef[l+1][j])+nf(l,j,ef[l][j]))*ef[l][j]
	    + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*dnf(l,j,ef[l][j])*(-ef[l][j+1]+ef[l][j])
	    + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*(nf(l,j+1,ef[l][j+1])+nf(l,j,ef[l][j]))*ef[l][j]
	    + lyr[k].grid2d[i][j-1].dSRH_Rec_ef(ef[l][j]);
	}
	else if(k<num_lyr-1 && i==lyr[k].numx()-1) {
	  printf("The interface is on the right\n");
	  if(chi(l,j)-chi(l+1,j)<-0.05) {
	    tmp1=-C[l][j]*(0.5*dx(l+1,j)+0.5*dx(l,j))*gamma*(vr/lyr[k].grid2d[i][jj].me())*(nf(l,j,ef[l][j])-nf(l,j,ef[l+1][j]));
	    tmp2=-C[l][j]*(0.5*dx(l+1,j)+0.5*dx(l,j))*gamma*(vr/lyr[k].grid2d[i][jj].me())*dnf(l,j,ef[l][j]);
	  }
	  else if(chi(l,j)-chi(l+1,j)>0.05){
	    tmp1= C[l][j]*(0.5*dx(l+1,j)+0.5*dx(l,j))*gamma*(vr/lyr[k+1].grid2d[0][jj].me())*(nf(l+1,j,ef[l+1][j])-nf(l+1,j,ef[l][j]));
	    tmp2=-C[l][j]*(0.5*dx(l+1,j)+0.5*dx(l,j))*gamma*(vr/lyr[k+1].grid2d[0][jj].me())*dnf(l+1,j,ef[l][j]);
	  }
	  else {
	    tmp1 = C[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*(nf(l+1,j,ef[l+1][j])+nf(l,j,ef[l][j]))*(-ef[l+1][j]+ef[l][j]);
	    tmp2 = C[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*dnf(l,j,ef[l][j])*(-ef[l+1][j]+ef[l][j])
	      + C[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*(nf(l+1,j,ef[l+1][j])+nf(l,j,ef[l][j]))*ef[l][j];
	  }
	  z = A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*(nf(l,j-1,ef[l][j-1])+nf(l,j,ef[l][j]))*(-ef[l][j-1]+ef[l][j])
	    + B[l][j]*0.5*(En[l-1][j]+En[l][j])*0.5*(nf(l-1,j,ef[l-1][j])+nf(l,j,ef[l][j]))*(-ef[l-1][j]+ef[l][j])
	    + tmp1
	    + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*(nf(l,j+1,ef[l][j+1])+nf(l,j,ef[l][j]))*(-ef[l][j+1]+ef[l][j])
	    + lyr[k].grid2d[i][jj].SRH_Rec_ef(ef[l][j]);
	  dz = A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*dnf(l,j,ef[l][j])*(-ef[l][j-1]+ef[l][j])
	    + A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*(nf(l,j-1,ef[l][j-1])+nf(l,j,ef[l][j]))*ef[l][j]
	    + B[l][j]*0.5*(En[l-1][j]+En[l][j])*0.5*dnf(l,j,ef[l][j])*(-ef[l-1][j]+ef[l][j])
	    + B[l][j]*0.5*(En[l-1][j]+En[l][j])*0.5*(nf(l-1,j,ef[l-1][j])+nf(l,j,ef[l][j]))*ef[l][j]
	    + tmp2
	    + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*dnf(l,j,ef[l][j])*(-ef[l][j+1]+ef[l][j])
	    + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*(nf(l,j+1,ef[l][j+1])+nf(l,j,ef[l][j]))*ef[l][j]
	    + lyr[k].grid2d[i][jj].dSRH_Rec_ef(ef[l][j]);
	}
	else if(k>0 && i==0) {
	  printf("The interface is on the left\n");
	  if(chi(l,j)-chi(l-1,j)>0.05) {
	    tmp1=-B[l][j]*(0.5*dx(l-1,j)+0.5*dx(l,j))*gamma*(vr/lyr[k-1].grid2d[lyr[k-1].numx()-1][jj].me())
	      *(nf(l-1,j,ef[l-1][j])-nf(l-1,j,ef[l][j]));
	    tmp2=B[l][j]*(0.5*dx(l-1,j)+0.5*dx(l,j))*gamma*(vr/lyr[k-1].grid2d[lyr[k-1].numx()-1][jj].me())*dnf(l-1,j,ef[l][j]);
	  }
	  else if(chi(l,j)-chi(l-1,j)<-0.05){
	    tmp1=B[l][j]*(0.5*dx(l-1,j)+0.5*dx(l,j))*gamma*(vr/lyr[k].grid2d[i][jj].me())*(nf(l,j,ef[l][j])-nf(l,j,ef[l-1][j]));
	    tmp2=B[l][j]*(0.5*dx(l-1,j)+0.5*dx(l,j))*gamma*(vr/lyr[k].grid2d[i][jj].me())*dnf(l,j,ef[l][j]);
	  }
	  else {
	    tmp1 = B[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*(nf(l+1,j,ef[l+1][j])+nf(l,j,ef[l][j]))*(-ef[l+1][j]+ef[l][j]);
	    tmp2 = B[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*dnf(l,j,ef[l][j])*(-ef[l+1][j]+ef[l][j])
	      + B[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*(nf(l+1,j,ef[l+1][j])+nf(l,j,ef[l][j]))*ef[l][j];
	  }
	  z = A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*(nf(l,j-1,ef[l][j-1])+nf(l,j,ef[l][j]))*(-ef[l][j-1]+ef[l][j])
	    + tmp1
	    + C[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*(nf(l+1,j,ef[l+1][j])+nf(l,j,ef[l][j]))*(-ef[l+1][j]+ef[l][j])
	    + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*(nf(l,j+1,ef[l][j+1])+nf(l,j,ef[l][j]))*(-ef[l][j+1]+ef[l][j])
	    + lyr[k].grid2d[i][j-1].SRH_Rec_ef(ef[l][j]);
	  dz = A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*dnf(l,j,ef[l][j])*(-ef[l][j-1]+ef[l][j])
	    + A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*(nf(l,j-1,ef[l][j-1])+nf(l,j,ef[l][j]))*ef[l][j]
	    + tmp2
	    + C[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*dnf(l,j,ef[l][j])*(-ef[l+1][j]+ef[l][j])
	    + C[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*(nf(l+1,j,ef[l+1][j])+nf(l,j,ef[l][j]))*ef[l][j]
	    + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*dnf(l,j,ef[l][j])*(-ef[l][j+1]+ef[l][j])
	    + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*(nf(l,j+1,ef[l][j+1])+nf(l,j,ef[l][j]))*ef[l][j]
	    + lyr[k].grid2d[i][j-1].dSRH_Rec_ef(ef[l][j]);
	}
	else if(k==0 && i==0) { /* Left contact */
	  printf("The left contact\n");
	  if(sim->barrier_left_contact!=0) {
	    l = 1;
	    z = A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*(lyr[k].grid2d[i][jj-1].nf(ef[l][j-1])+lyr[k].grid2d[i][jj].nf(ef[l][j]))
	      *(-ef[l][j-1]+ef[l][j])
	      + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*(lyr[k].grid2d[i][jj+1].nf(ef[l][j+1])+lyr[k].grid2d[i][jj].nf(ef[l][j]))
	      *(-ef[l][j+1]+ef[l][j])
	      + B[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*(lyr[k].grid2d[i+1][jj].nf(ef[l+1][j])+lyr[k].grid2d[i][jj].nf(ef[l][j]))
	      *(-ef[l+1][j]+ef[l][j])
	      + gamma*(vr/lyr[k].grid2d[i][jj].me())*(lyr[k].grid2d[i][jj].nf(ef[l][j])-lyr[k].grid2d[i][jj].nf(-sim->barrier_left_contact))
	      + lyr[k].grid2d[i][jj].SRH_Rec_ef(ef[l][j]);
	    dz = A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*lyr[k].grid2d[i][jj].dnf(ef[l][j])*(-ef[l][j-1]+ef[l][j])
	      + A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*(lyr[k].grid2d[i][jj-1].nf(ef[l][j-1])+lyr[k].grid2d[i][jj].nf(ef[l][j]))*ef[l][j]
	      + B[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*lyr[k].grid2d[i][jj].dnf(ef[l][j])*(-ef[l+1][j]+ef[l][j])
	      + B[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*(lyr[k].grid2d[i+1][jj].nf(ef[l+1][j])+lyr[k].grid2d[i][jj].nf(ef[l][j]))*ef[l][j]
	      + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*lyr[k].grid2d[i][jj].dnf(ef[l][j])*(-ef[l][j+1]+ef[l][j])
	      + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*(lyr[k].grid2d[i][jj+1].nf(ef[l][j+1])+lyr[k].grid2d[i][jj].nf(ef[l][j]))*ef[l][j]
	      + gamma*(vr/lyr[k].grid2d[i][jj].me())*lyr[k].grid2d[i][jj].dnf(ef[l][j])
	      + lyr[k].grid2d[i][jj].dSRH_Rec_ef(ef[l][j]);
	  }
	  else {
	    z = 0;
	    dz = 0;
	  }
	}
	else if(k==num_lyr-1 && i==lyr[k].numx()-1) { /* Right contact */
	  printf("The right contact\n");
	  if(sim->barrier_right_contact!=0) {
	    l = LayerToMesh(k,i);
	    z = A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*(lyr[k].grid2d[i][jj-1].nf(ef[l][j-1])+lyr[k].grid2d[i][jj].nf(ef[l][j]))
	      *(-ef[l][j-1]+ef[l][j])
	      + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*(lyr[k].grid2d[i][jj+1].nf(ef[l][j+1])+lyr[k].grid2d[i][jj].nf(ef[l][j]))
	      *(-ef[l][j+1]+ef[l][j])
	      + C[l][j]*0.5*(En[l-1][j]+En[l][j])*0.5*(lyr[k].grid2d[i-1][jj].nf(ef[l-1][j])+lyr[k].grid2d[i][jj].nf(ef[l][j]))
	      *(-ef[l-1][j]+ef[l][j])
	      - gamma*(vr/lyr[k].grid2d[i][jj].me())*(lyr[k].grid2d[i][jj].nf(ef[l][j])-lyr[k].grid2d[i][jj].nf(-sim->barrier_right_contact))
	      + lyr[k].grid2d[i][jj].SRH_Rec_ef(ef[l][j]);
	    dz = A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*lyr[k].grid2d[i][jj].dnf(ef[l][j])*(-ef[l][j-1]+ef[l][j])
	      + A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*(lyr[k].grid2d[i][jj-1].nf(ef[l][j-1])+lyr[k].grid2d[i][jj].nf(ef[l][j]))*ef[l][j]
	      + C[l][j]*0.5*(En[l-1][j]+En[l][j])*0.5*lyr[k].grid2d[i][jj].dnf(ef[l][j])*(-ef[l-1][j]+ef[l][j])
	      + C[l][j]*0.5*(En[l-1][j]+En[l][j])*0.5*(lyr[k].grid2d[i-1][jj].nf(ef[l-1][j])+lyr[k].grid2d[i][jj].nf(ef[l][j]))*ef[l][j]
	      + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*lyr[k].grid2d[i][jj].dnf(ef[l][j])*(-ef[l][j+1]+ef[l][j])
	      + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*(lyr[k].grid2d[i][jj+1].nf(ef[l][j+1])+lyr[k].grid2d[i][jj].nf(ef[l][j]))*ef[l][j]
	      - gamma*(vr/lyr[k].grid2d[i][jj].me())*lyr[k].grid2d[i][jj].dnf(ef[l][j])
	      + lyr[k].grid2d[i][jj].dSRH_Rec_ef(ef[l][j]);
	  }
	}
	else {
	  z =0;
	  dz = 0;
	}
	
	old = ef[l][j];
	if(z!=0 && fabs(dz)>1e-20) ef[l][j]=-z/dz;
	else if(z!=0 && fabs(dz)<=1e-20) {
	  fopen_s(&fp,"OUTPUT_FILES/e_zdependance.dat","w");
	  printf("# Point (%d,%d) in layer %d\n",l,j,k);
	  fprintf(fp,"# Point (%d,%d) in layer %d\n",i,j,k);
	  fprintf(fp,"# x  Y(x) and dY(x)/dx\n");
	  x = lyr[k].grid2d[i][jj].ret_phi_n()-5;
	  do {
	  z = A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*(nf(l,j-1,ef[l][j-1])+nf(l,j,x))*(-ef[l][j-1]+x)
	    + B[l][j]*0.5*(En[l-1][j]+En[l][j])*0.5*(nf(l-1,j,ef[l-1][j])+nf(l,j,x))*(-ef[l-1][j]+x)
	    + C[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*(nf(l+1,j,ef[l+1][j])+nf(l,j,x))*(-ef[l+1][j]+x)
	    + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*(nf(l,j+1,ef[l][j+1])+nf(l,j,x))*(-ef[l][j+1]+x)
	    + lyr[k].grid2d[i][j-1].SRH_Rec_ef(x);
	  dz = A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*dnf(l,j,x)*(-ef[l][j-1]+x)
	    + A[l][j]*0.5*(En[l][j-1]+En[l][j])*0.5*(nf(l,j-1,ef[l][j-1])+nf(l,j,x))*x
	    + B[l][j]*0.5*(En[l-1][j]+En[l][j])*0.5*dnf(l,j,x)*(-ef[l-1][j]+x)
	    + B[l][j]*0.5*(En[l-1][j]+En[l][j])*0.5*(nf(l-1,j,ef[l-1][j])+nf(l,j,x))*x
	    + C[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*dnf(l,j,x)*(-ef[l+1][j]+x)
	    + C[l][j]*0.5*(En[l+1][j]+En[l][j])*0.5*(nf(l+1,j,ef[l+1][j])+nf(l,j,x))*x
	    + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*dnf(l,j,x)*(-ef[l][j+1]+x)
	    + D[l][j]*0.5*(En[l][j+1]+En[l][j])*0.5*(nf(l,j+1,ef[l][j+1])+nf(l,j,x))*x
	    + lyr[k].grid2d[i][j-1].dSRH_Rec_ef(x);
	    fprintf(fp,"%lf  %le  %le\n",x,z,dz);
	    x +=0.01;
	  } while(x<lyr[k].grid2d[i][jj].ret_phi_n()+5);
	  fclose(fp);
	  exit(1);
	}
	iter++;
      } while(fabs(ef[l][j]-old)>PRECISION && iter<1000);
      if(iter>=1000) nrerror("2Dcontinuity: Precision not reached!");
    }
  }
  
  for(i=1;i<=xnum_mesh_points;i++) {
    ef[i][1]=ef[i][3];
    ef[i][ynum_mesh_points]=ef[i][ynum_mesh_points-2];
  }
  
  printf("Holes...\n");
  
  /* HOLE SOLVER */
  
  for(k=0;k<num_lyr;k++) {
    
    num_x = lyr[k].numx();
    num_y = lyr[k].numy();
    
    for(j=1;j<num_y-1;j++)
      for(i=0;i<num_x;i++) {
	l = LayerToMesh(k,i)-1;
	iter = 0;
	do {
	  if((k>0 && k<num_lyr-1 && i>0 && i<num_x-1) 
	     || (k==0 && i>0 && i<num_x-1 ) 
	     || (k==num_lyr-1 && i>0 && i<num_x-1 ) ) {
	    z = A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j-1].pf(hf[l][j-1])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l][j-1]+hf[l][j])
	      + B[l][j]*0.5*(Ep[l-1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i-1][j].pf(hf[l-1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l-1][j]+hf[l][j])
	      + C[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i+1][j].pf(hf[l+1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l+1][j]+hf[l][j])
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j+1].pf(hf[l][j+1])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l][j+1]+hf[l][j])
	      + lyr[k].grid2d[i][j].SRH_Rec_hf(hf[i][j]);
	    dz = A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l][j-1]+hf[l][j])
	      + A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j-1].pf(hf[l][j-1])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + B[l][j]*0.5*(Ep[l-1][j]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l-1][j]+hf[l][j])
	      + B[l][j]*0.5*(Ep[l-1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i-1][j].pf(hf[l-1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + C[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l+1][j]+hf[l][j])
	      + C[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i+1][j].pf(hf[l+1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l][j+1]+hf[l][j])
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j+1].pf(hf[l][j+1])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + lyr[k].grid2d[i][j].dSRH_Rec_hf(hf[i][j]);
	  }
	  else if(k<num_lyr-1 && i==num_x-1) {
	    if(lyr[k].grid2d[i][j].chi()-lyr[k+1].grid2d[0][j].chi()<-0.05) {
	      tmp1=-C[l][j]*(0.5*lyr[k+1].grid2d[0][j].xsize()+0.5*lyr[k].grid2d[i][j].xsize())
		*gamma*(vr/lyr[k].grid2d[i][j].me())*(lyr[k].grid2d[i][j].pf(hf[l][j])-lyr[k].grid2d[i][j].pf(hf[l+1][j]));
	      tmp2=-C[l][j]*(0.5*lyr[k+1].grid2d[0][j].xsize()+0.5*lyr[k].grid2d[i][j].xsize())
		*gamma*(vr/lyr[k].grid2d[i][j].me())*lyr[k].grid2d[i][j].dpf(hf[l][j]);
	    }
	    else if(lyr[k].grid2d[i][j].chi()-lyr[k+1].grid2d[0][j].chi()>0.05){
	      tmp1=C[l][j]*(0.5*lyr[k+1].grid2d[0][j].xsize()+0.5*lyr[k].grid2d[i][j].xsize())
		*gamma*(vr/lyr[k+1].grid2d[0][j].me())*(lyr[k+1].grid2d[0][j].pf(hf[l+1][j])-lyr[k+1].grid2d[0][j].pf(hf[l][j]));
	      tmp2=-C[l][j]*(0.5*lyr[k+1].grid2d[0][j].xsize()+0.5*lyr[k].grid2d[i][j].xsize())
		*gamma*(vr/lyr[k+1].grid2d[0][j].me())*lyr[k+1].grid2d[0][j].dpf(hf[l][j]);
	    }
	    else {
	      tmp1 = C[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*(lyr[k+1].grid2d[0][j].pf(hf[l+1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))
		*(-hf[l+1][j]+hf[l][j]);
	      tmp2 = C[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l+1][j]+hf[l][j])
		+ C[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*(lyr[k+1].grid2d[0][j].pf(hf[l+1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j];
	    }

	    z = A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j-1].pf(hf[l][j-1])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l][j-1]+hf[l][j])
	      + B[l][j]*0.5*(Ep[l-1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i-1][j].pf(hf[l-1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l-1][j]+hf[l][j])
	      + tmp1
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j+1].pf(hf[l][j+1])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l][j+1]+hf[l][j])
	      + lyr[k].grid2d[i][j].SRH_Rec_hf(hf[i][j]);
	    dz = A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l][j-1]+hf[l][j])
	      + A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j-1].pf(hf[l][j-1])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + B[l][j]*0.5*(Ep[l-1][j]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l-1][j]+hf[l][j])
	      + B[l][j]*0.5*(Ep[l-1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i-1][j].pf(hf[l-1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + tmp2
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l][j+1]+hf[l][j])
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j+1].pf(hf[l][j+1])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + lyr[k].grid2d[i][j].dSRH_Rec_hf(hf[i][j]);

	  }
	  else if(k>0 && i==0) {

	    if(lyr[k].grid2d[i][j].chi()-lyr[k-1].grid2d[lyr[k-1].numx()-1][j].chi()>0.05) {
	      tmp1=-B[l][j]*(0.5*lyr[k-1].grid2d[lyr[k-1].numx()-1][j].xsize()+0.5*lyr[k].grid2d[i][j].xsize())
		*gamma*(vr/lyr[k-1].grid2d[lyr[k-1].numx()-1][j].me())
		*(lyr[k-1].grid2d[lyr[k-1].numx()-1][j].pf(hf[l-1][j])-lyr[k-1].grid2d[lyr[k-1].numx()-1][j].pf(hf[l][j]));
	      tmp2=B[l][j]*(0.5*lyr[k-1].grid2d[lyr[k-1].numx()-1][j].xsize()+0.5*lyr[k].grid2d[i][j].xsize())
		*gamma*(vr/lyr[k-1].grid2d[lyr[k-1].numx()-1][j].me())*lyr[k-1].grid2d[lyr[k-1].numx()-1][j].dpf(hf[l][j]);
	    }
	    else if(lyr[k].grid2d[i][j].chi()-lyr[k-1].grid2d[lyr[k-1].numx()-1][j].chi()<-0.05){
	      tmp1=B[l][j]*(0.5*lyr[k-1].grid2d[lyr[k-1].numx()-1][j].xsize()+0.5*lyr[k].grid2d[i][j].xsize())
		*gamma*(vr/lyr[k].grid2d[i][j].me())*(lyr[k].grid2d[i][j].pf(hf[l][j])-lyr[k].grid2d[i][j].pf(hf[l-1][j]));
	      tmp2=B[l][j]*(0.5*lyr[k-1].grid2d[lyr[k-1].numx()-1][j].xsize()+0.5*lyr[k].grid2d[i][j].xsize())
		*gamma*(vr/lyr[k].grid2d[i][j].me())*lyr[k].grid2d[i][j].dpf(hf[l][j]);
	    }
	    else {
	      tmp1 = B[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*(lyr[k+1].grid2d[0][j].pf(hf[l+1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))
		*(-hf[l+1][j]+hf[l][j]);
	      tmp2 = B[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l+1][j]+hf[l][j])
		+ B[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*(lyr[k+1].grid2d[0][j].pf(hf[l+1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j];
	    }
	    z = A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j-1].pf(hf[l][j-1])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l][j-1]+hf[l][j])
	      + tmp1
	      + C[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i+1][j].pf(hf[l+1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l+1][j]+hf[l][j])
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j+1].pf(hf[l][j+1])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l][j+1]+hf[l][j])
	      + lyr[k].grid2d[i][j].SRH_Rec_hf(hf[i][j]);
	    dz = A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l][j-1]+hf[l][j])
	      + A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j-1].pf(hf[l][j-1])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + tmp2
	      + C[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l+1][j]+hf[l][j])
	      + C[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i+1][j].pf(hf[l+1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l][j+1]+hf[l][j])
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j+1].pf(hf[l][j+1])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + lyr[k].grid2d[i][j].dSRH_Rec_hf(hf[i][j]);
	
	  }
	  else if(k==0 && i==0) { /* Left contact */
	    if(sim->barrier_left_contact!=0) {
	      l = 1;
	      z = A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j-1].pf(hf[l][j-1])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l][j-1]+hf[l][j])
		+ D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j+1].pf(hf[l][j+1])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l][j+1]+hf[l][j])
		+ B[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i+1][j].pf(hf[l+1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l+1][j]+hf[l][j])
		+ gamma*(vr/lyr[k].grid2d[i][j].me())*(lyr[k].grid2d[i][j].pf(hf[l][j])-lyr[k].grid2d[i][j].pf(-sim->barrier_left_contact))
		+ lyr[k].grid2d[i][j].SRH_Rec_hf(hf[i][j]);
	      dz = A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l][j-1]+hf[l][j])
	      + A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j-1].pf(hf[l][j-1])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + B[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l+1][j]+hf[l][j])
	      + B[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i+1][j].pf(hf[l+1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l][j+1]+hf[l][j])
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j+1].pf(hf[l][j+1])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + gamma*(vr/lyr[k].grid2d[i][j].me())*lyr[k].grid2d[i][j].dpf(hf[l][j])
	      + lyr[k].grid2d[i][j].dSRH_Rec_hf(hf[i][j]);
	    }
	    else {
	      z = 0;
	      dz = 0;
	    }
	  }
	  else if(k==num_lyr-1 && i==num_x-1) { /* Right contact */
	    if(sim->barrier_right_contact!=0) {
	      l = LayerToMesh(k,i);
	      z = A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j-1].pf(hf[l][j-1])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l][j-1]+hf[l][j])
		+ D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j+1].pf(hf[l][j+1])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l][j+1]+hf[l][j])
		+ C[l][j]*0.5*(Ep[l-1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i-1][j].pf(hf[l-1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))
	      *(-hf[l-1][j]+hf[l][j])
		- gamma*(vr/lyr[k].grid2d[i][j].me())*(lyr[k].grid2d[i][j].pf(hf[l][j])-lyr[k].grid2d[i][j].pf(-sim->barrier_right_contact))
		+ lyr[k].grid2d[i][j].SRH_Rec_hf(hf[i][j]);
	      dz = A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l][j-1]+hf[l][j])
	      + A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j-1].pf(hf[l][j-1])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + C[l][j]*0.5*(Ep[l-1][j]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l-1][j]+hf[l][j])
	      + C[l][j]*0.5*(Ep[l-1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i-1][j].pf(hf[l-1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l][j+1]+hf[l][j])
	      + D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j+1].pf(hf[l][j+1])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
	      - gamma*(vr/lyr[k].grid2d[i][j].me())*lyr[k].grid2d[i][j].dpf(hf[l][j])
	      + lyr[k].grid2d[i][j].dSRH_Rec_hf(hf[i][j]);
	    }
	  }
	  else {
	    z =0;
	    dz = 0;
	  }
	  
	  old = hf[i][j];
	  if(z!=0 && fabs(dz)>1e-20) hf[i][j]=-z/dz;
	  else if(z!=0 && fabs(dz)<=1e-20) {
	    fopen_s(&fp,"OUTPUT_FILES/h_zdependance.dat","w");
	    fprintf(fp,"# Point (%d,%d) in layer %d\n",i,j,k);
	    fprintf(fp,"# x  Y(x) and dY(x)/dx\n");
	    x = lyr[k].grid2d[i][j].ret_phi_p()-5;
	    do {
	      z = A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j-1].pf(hf[l][j-1])+lyr[k].grid2d[i][j].pf(hf[l][j]))
		*(-hf[l][j-1]+hf[l][j])
		+ B[l][j]*0.5*(Ep[l-1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i-1][j].pf(hf[l-1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))
		*(-hf[l-1][j]+hf[l][j])
		+ C[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i+1][j].pf(hf[l+1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))
		*(-hf[l+1][j]+hf[l][j])
		+ D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j+1].pf(hf[l][j+1])+lyr[k].grid2d[i][j].pf(hf[l][j]))
		*(-hf[l][j+1]+hf[l][j])
		+ lyr[k].grid2d[i][j].SRH_Rec_hf(hf[i][j]);
	      dz = A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l][j-1]+hf[l][j])
		+ A[l][j]*0.5*(Ep[l][j-1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j-1].pf(hf[l][j-1])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
		+ B[l][j]*0.5*(Ep[l-1][j]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l-1][j]+hf[l][j])
		+ B[l][j]*0.5*(Ep[l-1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i-1][j].pf(hf[l-1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
		+ C[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l+1][j]+hf[l][j])
		+ C[l][j]*0.5*(Ep[l+1][j]+Ep[l][j])*0.5*(lyr[k].grid2d[i+1][j].pf(hf[l+1][j])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
		+ D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*lyr[k].grid2d[i][j].dpf(hf[l][j])*(-hf[l][j+1]+hf[l][j])
		+ D[l][j]*0.5*(Ep[l][j+1]+Ep[l][j])*0.5*(lyr[k].grid2d[i][j+1].pf(hf[l][j+1])+lyr[k].grid2d[i][j].pf(hf[l][j]))*hf[l][j]
		+ lyr[k].grid2d[i][j].dSRH_Rec_hf(hf[i][j]);
	      fprintf(fp,"%lf  %le  %le",x,z,dz);
	      x +=0.01;
	    } while(x<lyr[k].grid2d[i][j].ret_phi_p()+5);
	    fclose(fp);
	  }
	  iter++;
	} while(fabs(hf[i][j]-old)>PRECISION && iter<1000);
	if(iter>=1000) nrerror("2Dcontinuity hole: Precision not reached!");
      }

    for(i=1;i<=num_x;i++) {
      hf[i][0]=hf[i][2];
      hf[i][num_y-1]=hf[i][num_y-1];
    }
    
  }

  for(j=0;j<ynum_mesh_points;j++)
    for(i=0;i<xnum_mesh_points;i++) {
      set_phin(i+1,j+1,ef[i][j]);
      set_phip(i+1,j+1,hf[i][j]);
    }

  free_matrix(A,2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  free_matrix(B,2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  free_matrix(C,2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  free_matrix(D,2,xnum_mesh_points-1,2,ynum_mesh_points-1);
  free_matrix(En,1,xnum_mesh_points,1,ynum_mesh_points);
  free_matrix(Ep,1,xnum_mesh_points,1,ynum_mesh_points);
  
  free_matrix(ef,1,xnum_mesh_points,1,ynum_mesh_points);
  free_matrix(hf,1,xnum_mesh_points,1,ynum_mesh_points);

}


/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */

double ** Device::ElectronDiscretization1D(short type, int & num) {

  int i;
  double a1,a2,a3;
  double ** Matrix;
  
  /* ALLOCATE MATRIX */

  num = xnum_mesh_points-2;
  Matrix = matrix(1,num,1,4);

  for(i=2;i<=num+1;i++) {

    /* left current coefficient */
    a1 = (phi(i)-phi(i-1))/VT + (chi(i)-chi(i-1))/VT + log(Nc(i)/Nc(i-1));
    a2 = Bernoulli(a1);
    a3 = 2.0*a2*Dn(i)*Nc(i)*exp((phi(i)+chi(i))/VT);
    Matrix[i-1][1]= 2.0*a3/(dx(i)+dx(i-1));

    /* right current coefficient */
    a1 = (phi(i)- phi(i+1))/VT + (chi(i)-chi(i+1))/VT + log(Nc(i)/Nc(i+1));
    a2  = Bernoulli(a1);
    a3	= 2.0*a2*Dn(i)*Nc(i)*exp((phi(i)+chi(i))/VT);
    Matrix[i-1][3]= 2.0*a3/(dx(i)+dx(i+1));

    /* middle coefficient */
    Matrix[i-1][2] =- Matrix[i-1][1] - Matrix[i-1][3];

    /* known coefficient */
    if(sim->traps) {
      Matrix[i-1][4] = ( (n(i)*p(i)-ni(i)*ni(i))/ ( Tn(i)*(p(i)+ni(i))+ Tp(i)*(n(i)+ni(i))) )*(2.0*dx(i)+dx(i-1)+dx(i+1))
	              +( (n(i)*p(i)-ni(i)*ni(i))/ ( Tntrap(i)*(p(i)+p1(i))+ Tptrap(i)*(n(i)+n1(i))) )*(2.0*dx(i)+dx(i-1)+dx(i+1));
    }
    else  Matrix[i-1][4] = ( (n(i)*p(i)-ni(i)*ni(i))/ ( Tn(i)*(p(i)+ni(i))+ Tp(i)*(n(i)+ni(i))) )*(2.0*dx(i)+dx(i-1)+dx(i+1));
  }

  /* LEFT BOUNDARY */

  Matrix[1][4] -= Matrix[1][1]*exp(-phin(1)/VT);
  Matrix[1][1] = 0;

  /* RIGHT BOUNDARY */

  Matrix[num][4] -= Matrix[num][3]*exp(-phin(num+2)/VT);
  Matrix[num][3] = 0;

  return Matrix;
}
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */

double ** Device::ElectronDiscretization1D_3(short type, int & num) {

  int i,k,j,j_tmp;
  double a1,a2,a3;
  double ** Matrix;
  const double STEP_BARRIER = 0.2; // (eV)
  const double MAX_QW_WIDTH = 5e-6; // (cm)
  const double COEFF_TH = 1e-4*(Q/PLANKUT)*(Q/PLANKUT)*(Q/PLANKUT)*VT*VT/(2.0*PI*PI);

  /* ALLOCATE MATRIX */

  num = xnum_mesh_points-2;
  Matrix = matrix(1,num,1,4);

  for(j=2;j<=num+1;j++) {
    MeshToLayer(j,k,i);
    if(k==0 && i==1) { /* Left Contact */
      if(type==1) { /* Ohmic */
	/* left current coefficient */
	a1 = (phi(2)-phi(1))/VT + (chi(2)-chi(1))/VT + log(Nc(2)/Nc(1));
	a2 = Bernoulli(a1);
	a3 = 2.0*a2*Dn(2)*Nc(2)*exp((phi(2)+chi(2))/VT);
	Matrix[1][1]= 2.0*a3/(dx(2)+dx(1));
	/* right current coefficient */
	a1 = (phi(2)- phi(3))/VT + (chi(2)-chi(3))/VT + log(Nc(2)/Nc(3));
	a2 = Bernoulli(a1);
	a3 = 2.0*a2*Dn(2)*Nc(2)*exp((phi(2)+chi(2))/VT);
	Matrix[1][3]= 2.0*a3/(dx(2)+dx(3));
	/* middle coefficient */
	Matrix[1][2] =- Matrix[1][1] - Matrix[1][3];
	/* known coefficient */
	if(sim->traps) Matrix[1][4] = ( SHR_Rec(2) )*(2.0*dx(2)+dx(1)+dx(3))+( Trap_Rec_e(2)+Trap_Rec_h(2) )*(2.0*dx(2)+dx(1)+dx(3));
	else  Matrix[1][4] = ( SHR_Rec(2) )*(2.0*dx(2)+dx(1)+dx(3));
	Matrix[1][4] -= Matrix[1][1]*exp(-phin(1)/VT);
	Matrix[1][1] = 0;
      }
      else { /* Schottky */
	printf("LEFT SCHOTTKY: NOT YET IMPLEMENTED\n");
	exit(1);
      }
    }
    else if(k==num_lyr-1 && i==lyr[k].numx()-2) { /* Right Contact */
      if(type==1) { /* Ohmic */
	/* left current coefficient */
	a1 = (phi(j)-phi(j-1))/VT + (chi(j)-chi(j-1))/VT + log(Nc(j)/Nc(j-1));
	a2 = Bernoulli(a1);
	a3 = 2.0*a2*Dn(j)*Nc(j)*exp((phi(j)+chi(j))/VT);
	Matrix[j-1][1]= 2.0*a3/(dx(j)+dx(j-1));
	
	/* right current coefficient */
	a1 = (phi(j)- phi(j+1))/VT + (chi(j)-chi(j+1))/VT + log(Nc(j)/Nc(j+1));
	a2 = Bernoulli(a1);
	a3 = 2.0*a2*Dn(j)*Nc(j)*exp((phi(j)+chi(j))/VT);
	Matrix[j-1][3]= 2.0*a3/(dx(j)+dx(j+1));
	
	/* middle coefficient */
	Matrix[j-1][2] =- Matrix[j-1][1] - Matrix[j-1][3];
	
	/* known coefficient */
	if(sim->traps) Matrix[j-1][4] = ( SHR_Rec(j)+Trap_Rec_e(j)+Trap_Rec_h(j)+Opt_Rec(j) )*(2.0*dx(j)+dx(j-1)+dx(j+1));
	else  Matrix[j-1][4] = SHR_Rec(j)*(2.0*dx(j)+dx(j-1)+dx(j+1));
	Matrix[j-1][4] -=  Matrix[j-1][3]*exp(-phin(j+1)/VT);
	Matrix[j-1][3] = 0.0;
      }
      else { /* Schottky */
	printf("RIGHT SCHOTTKY: NOT YET IMPLEMENTED\n");
	exit(1);
      }
    }
    else if(fabs(chi(j)-chi(j+1))>STEP_BARRIER) { /* Step on the right! */
      
      if(chi(j)-chi(j+1)>STEP_BARRIER) { // STEP-UP
	if(k>0 && ((lyr[k].xsize()<MAX_QW_WIDTH && (lyr[k].grid1d[0].chi()-lyr[k-1].grid1d[lyr[k-1].numx()-1].chi())>STEP_BARRIER) 
		   || lyr[k].grid1d[i].qptr_e!=0)) { // OUT-OF-A-QUANTUM-WELL
	  /* left current coefficient */
	  a1 = (phi(j)-phi(j-1))/VT;
	  a2 = Bernoulli(a1);
	  a3 = 2.0*a2*Dn(j)*Nc(j+1)*exp((phi(j)+chi(j+1))/VT);
	  Matrix[j-1][1]= 2.0*a3/(dx(j)+dx(j-1));
	  
	  /* right current coefficient */
	  a1 = (phi(j)- phi(j+1))/VT;
	  a2 = Bernoulli(a1);
	  a3 = 2.0*a2*Dn(j)*Nc(j+1)*exp((phi(j)+chi(j+1))/VT);
	  Matrix[j-1][3]= 2.0*a3/(dx(j)+dx(j+1));
	  
	  /* middle coefficient */
	  Matrix[j-1][2] =- Matrix[j-1][1] - Matrix[j-1][3];
	  
	  /* known coefficient */
	  if(sim->traps) Matrix[j-1][4] = (QW_Rec(j) + SHR_Rec(j)+Trap_Rec_e(j)+Trap_Rec_h(j)+Opt_Rec(j+1) )*(2.0*dx(j)+dx(j-1)+dx(j+1));
	  else  Matrix[j-1][4] = (QW_Rec(j) + SHR_Rec(j))*(2.0*dx(j)+dx(j-1)+dx(j+1));	 
	  
	}
	else { // JUST-A-BARRIER
	  /* left current coefficient */
	  a1 = (phi(j)-phi(j-1))/VT+ (chi(j)-chi(j-1))/VT + log(Nc(j)/Nc(j-1));
	  a2 = Bernoulli(a1);
	  a3 = 2.0*a2*Dn(j)*Nc(j+1)*exp((phi(j)+chi(j+1))/VT);
	  Matrix[j-1][1]= 2.0*a3/(dx(j)+dx(j-1));
	  
	  /* right current coefficient */
	  a3 = 2.0/((1.0/me(j))+(1.0/me(j+1)));
	  Matrix[j-1][3]=COEFF_TH*a3*exp((phi(j)+chi(j+1))/VT);
	  
	  /* middle coefficient */
	  Matrix[j-1][2] =- Matrix[j-1][1] - Matrix[j-1][3];
	  
	  /* known coefficient */
	  if(sim->traps) Matrix[j-1][4] = (SHR_Rec(j)+Trap_Rec_e(j)+Trap_Rec_h(j)+Opt_Rec(j+1) )*(2.0*dx(j)+dx(j-1)+dx(j+1));
	  else  Matrix[j-1][4] = SHR_Rec(j)*(2.0*dx(j)+dx(j-1)+dx(j+1));
	}
	
      }
      else { // STEP-DOWN
	if(k<num_lyr-1 && ((lyr[k].xsize()<MAX_QW_WIDTH && (lyr[k].grid1d[lyr[k].numx()-1].chi()-lyr[k+1].grid1d[0].chi())>STEP_BARRIER) 
			   || lyr[k].grid1d[i].qptr_e!=0)) { // IN-A-QUANTUM-WELL
	  
	  /* left current coefficient */
	  a1 = (phi(j)-phi(j-1))/VT;
	  a2 = Bernoulli(a1);
	  a3 = 2.0*a2*Dn(j)*Nc(j-1)*exp((phi(j)+chi(j-1))/VT);
	  Matrix[j-1][1]= 2.0*a3/(dx(j)+dx(j-1));
	  
	  /* right current coefficient */
	  a1 = (phi(j)- phi(j+1))/VT;
	  a2 = Bernoulli(a1);
	  a3 = 2.0*a2*Dn(j)*Nc(j-1)*exp((phi(j)+chi(j-1))/VT);
	  Matrix[j-1][3]= 2.0*a3/(dx(j)+dx(j+1));
	  
	  /* middle coefficient */
	  Matrix[j-1][2] =- Matrix[j-1][1] - Matrix[j-1][3];
	  
	  /* known coefficient */
	  if(sim->traps) Matrix[j-1][4] = (QW_Rec(j) + SHR_Rec(j)+Trap_Rec_e(j)+Trap_Rec_h(j)+Opt_Rec(j-1) )*(2.0*dx(j)+dx(j-1)+dx(j+1));
	  else  Matrix[j-1][4] = (QW_Rec(j) + SHR_Rec(j))*(2.0*dx(j)+dx(j-1)+dx(j+1));	 
	  
	}
	else { // JUST-A-JUMP-DOWN
	  /* left current coefficient */
	  a1 = (phi(j)-phi(j-1))/VT+ (chi(j)-chi(j-1))/VT + log(Nc(j)/Nc(j-1));
	  a2 = Bernoulli(a1);
	  a3 = 2.0*a2*Dn(j)*Nc(j+1)*exp((phi(j)+chi(j+1))/VT);
	  Matrix[j-1][1]= 2.0*a3/(dx(j)+dx(j-1));
	  
	  /* right current coefficient */
	  a3 = 2.0/((1.0/me(j))+(1.0/me(j+1)));
	  Matrix[j-1][3]=COEFF_TH*a3*exp((phi(j)+chi(j))/VT);
	  
	  /* middle coefficient */
	  Matrix[j-1][2] =- Matrix[j-1][1] - Matrix[j-1][3];
	  
	  /* known coefficient */
	  if(sim->traps) Matrix[j-1][4] = (SHR_Rec(j)+Trap_Rec_e(j)+Trap_Rec_h(j)+Opt_Rec(j+1) )*(2.0*dx(j)+dx(j-1)+dx(j+1));
	  else  Matrix[j-1][4] = SHR_Rec(j)*(2.0*dx(j)+dx(j-1)+dx(j+1));	
	}
	}
    }
    else {

      if( (lyr[k].xsize()<MAX_QW_WIDTH) && ((lyr[k].grid1d[lyr[k].numx()-1].chi()-lyr[k+1].grid1d[0].chi())>STEP_BARRIER) &&
	  ((lyr[k].grid1d[0].chi()-lyr[k-1].grid1d[lyr[k-1].numx()-1].chi())>STEP_BARRIER) ) { // IN THE QW
	j_tmp = LayerToMesh(k-1,lyr[k-1].numx()-1);
	/* left current coefficient */
	a1 = (phi(j)-phi(j-1))/VT;
	a2 = Bernoulli(a1);
	a3 = 2.0*a2*Dn(j)*Nc(j_tmp)*exp((phi(j)+chi(j_tmp))/VT);
	Matrix[j-1][1]= 2.0*a3/(dx(j)+dx(j-1));
	
	/* right current coefficient */
	a1 = (phi(j)- phi(j+1))/VT;
	a2 = Bernoulli(a1);
	a3 = 2.0*a2*Dn(j)*Nc(j_tmp)*exp((phi(j)+chi(j_tmp))/VT);
	Matrix[j-1][3]= 2.0*a3/(dx(j)+dx(j+1));
	
	/* middle coefficient */
	Matrix[j-1][2] =- Matrix[j-1][1] - Matrix[j-1][3];
	
	/* known coefficient */
	if(sim->traps) Matrix[j-1][4] = ( QW_Rec(j)+SHR_Rec(j)+Trap_Rec_e(j)+Trap_Rec_h(j)+Opt_Rec(j) )*(2.0*dx(j)+dx(j-1)+dx(j+1));
	else  Matrix[j-1][4] = (QW_Rec(j)+SHR_Rec(j))*(2.0*dx(j)+dx(j-1)+dx(j+1));	
      }
      else { // IN-THE-BULK
	/* left current coefficient */
	a1 = (phi(j)-phi(j-1))/VT + (chi(j)-chi(j-1))/VT + log(Nc(j)/Nc(j-1));
	a2 = Bernoulli(a1);
	a3 = 2.0*a2*Dn(j)*Nc(j)*exp((phi(j)+chi(j))/VT);
	Matrix[j-1][1]= 2.0*a3/(dx(j)+dx(j-1));
	
	/* right current coefficient */
	a1 = (phi(j)- phi(j+1))/VT + (chi(j)-chi(j+1))/VT + log(Nc(j)/Nc(j+1));
	a2 = Bernoulli(a1);
	a3 = 2.0*a2*Dn(j)*Nc(j)*exp((phi(j)+chi(j))/VT);
	Matrix[j-1][3]= 2.0*a3/(dx(j)+dx(j+1));
	
	/* middle coefficient */
	Matrix[j-1][2] =- Matrix[j-1][1] - Matrix[j-1][3];
	
	/* known coefficient */
	if(sim->traps) Matrix[j-1][4] = ( SHR_Rec(j)+Trap_Rec_e(j)+Trap_Rec_h(j)+Opt_Rec(j) )*(2.0*dx(j)+dx(j-1)+dx(j+1));
	else  Matrix[j-1][4] = SHR_Rec(j)*(2.0*dx(j)+dx(j-1)+dx(j+1));	
      }
    }
  }

  return Matrix;
}

/* --------------------------------------------------------------- */

 double ** Device::HoleDiscretization1D(short type, int & num) {

  int i;
  double a1,a2,a3;
  double ** Matrix;
  

  /* ALLOCATE MATRIX */

  num = xnum_mesh_points-2;
  Matrix = matrix(1,num,1,4);

  for(i=2;i<=num+1;i++) {

    /* left current coefficient */
    a1 = (-phi(i)+phi(i-1))/VT + (-chi(i)+chi(i-1))/VT + (-Eg(i)+Eg(i-1))/VT + log(Nv(i)/Nv(i-1));
    a2 = Bernoulli(a1);
    a3 = 2.0*a2*Dp(i)*Nv(i)*exp(-(phi(i)+chi(i)+Eg(i))/VT);
    Matrix[i-1][1]= 2.0*a3/(dx(i)+dx(i-1));

    /* right current coefficient */
    a1 = (-phi(i)+phi(i+1))/VT + (-chi(i)+chi(i+1))/VT + (-Eg(i)+Eg(i+1))/VT + log(Nv(i)/Nv(i+1));
    a2 = Bernoulli(a1);
    a3 = 2.0*a2*Dp(i)*Nv(i)*exp(-(phi(i)+chi(i)+Eg(i))/VT);
    Matrix[i-1][3]= 2.0*a3/(dx(i)+dx(i+1));

    /* middle coefficient */
    Matrix[i-1][2] = - Matrix[i-1][1] - Matrix[i-1][3];

    /* known coefficient */
    if(sim->traps) {
      Matrix[i-1][4] = ( (n(i)*p(i)-ni(i)*ni(i))/ ( Tn(i)*(p(i)+ni(i))+ Tp(i)*(n(i)+ni(i))) )*(2.0*dx(i)+dx(i-1)+dx(i+1))
	              +( (n(i)*p(i)-ni(i)*ni(i))/ ( Tntrap(i)*(p(i)+p1(i))+ Tptrap(i)*(n(i)+n1(i))) )*(2.0*dx(i)+dx(i-1)+dx(i+1));
    }
    else  Matrix[i-1][4] = ( (n(i)*p(i)-ni(i)*ni(i))/ ( Tn(i)*(p(i)+ni(i))+ Tp(i)*(n(i)+ni(i))) )*(2.0*dx(i)+dx(i-1)+dx(i+1));

  }

  /* LEFT BOUNDARY */

  Matrix[1][4] -= Matrix[1][1]*exp(phip(1)/VT);
  Matrix[1][1] = 0;    

  /* RIGHT BOUNDARY */

  Matrix[num][4] -= Matrix[num][3]*exp(phip(num+2)/VT);
  Matrix[num][3] = 0;

  return Matrix;
}

/* ------------------------------------------------------------ */

void Device::ShowMatrix(double ** Matrix, int nrow, int ncol) {
  int i,j;
  FILE *fp;

  fopen_s(&fp,"OUTPUT_FILES/matrix.dat","w");
  if(ncol>nrow) {
    printf("I have transposed the matrix\n");
    for(i=1;i<=ncol;i++) {
      for(j=1;j<=nrow;j++) fprintf(fp,"%le  ",Matrix[j][i]);
      fprintf(fp,"\n");
    }
  }
  else {
    for(i=1;i<=nrow;i++) {
      for(j=1;j<=ncol;j++) fprintf(fp,"%le  ",Matrix[i][j]);
      fprintf(fp,"\n");
    }
  }
  fclose(fp);
  printf("File matrix.dat dumped\n");
  
}

/* ------------------------------------------------------------ */

void Device::ShowVector(double * Vector, int n) {
  int j;
  FILE *fp;

  fopen_s(&fp,"OUTPUT_FILES/vector.dat","w");
  for(j=1;j<=n;j++) fprintf(fp,"%le\n",Vector[j]);
  fclose(fp);
  printf("File vector.dat dumped\n");
  
}
/* ------------------------------------------------------------ */

 double Device::ElectronSolver(short type, double x, short debug) {

  int n,j,i;
  double error,tmp;
  double ** Matrix = ElectronDiscretization1D(1,n);
  double * u = vector(1,n);

  if(debug) printf("Electron Solver...\n");

  if(type==0) tridag(Matrix[1],Matrix[2],Matrix[3],Matrix[4],u,n);
  else if(type==1) {
    double ** a = matrix(1,n,1,3);
    double ** al = matrix(1,n,1,1);
    int * indx = ivector(1,n);
    double d;
    for(i=1;i<=n;i++) for(j=1;j<=3;j++) a[i][j] = Matrix[i][j];
    for(i=1;i<=n;i++) u[i] = Matrix[i][4];
    bandec(a,n,1,1,al,indx,&d);
    banbks(a,n,1,1,al,indx,u);
    free_ivector(indx,1,n);
    free_matrix(a,1,n,1,3);
    free_matrix(al,1,n,1,1);
  }
  else if (type==2) {
    for(i=1;i<=n;i++) u[i] = exp(-phin(i+1)/VT);
    u[1] = (Matrix[1][4]-Matrix[1][3]*u[2])/Matrix[1][2];
    for(i=2;i<n;i++) u[i]= (Matrix[i][4] -Matrix[i][1]*u[i-1]-Matrix[i][3]*u[i+1])/Matrix[i][2];
    u[n] = (Matrix[n][4]-Matrix[n][1]*u[n-1])/Matrix[n][2];
  }
  else nrerror("Select right electron solver type!");
  
  // if(debug) ShowVector(u,n);

  error = 0.0;
  for(i=1;i<=n;i++) {
    tmp = fabs(phin(i+1)+VT*log(u[i]));
    if(tmp>1e3) {
      printf("%d: phin=%le  new=%le\n",i,phin(i+1),-VT*log(u[i]));
    }
    error += tmp;
    set_phin(i+1,(1.0-x)*phin(i+1) + x*(-VT*log(u[i])));     
  }

  free_matrix(Matrix,1,n,1,4);
  free_vector(u,1,n);

  if(debug) printf("Electron solver done!\n");

  return error;

}

/* ------------------------------------------------------------ */

 double Device::HoleSolver(short type, double x, short debug) {

  int n,j,i;
  double error,tmp;
  double ** Matrix = HoleDiscretization1D(0,n);
  double * v = vector(1,n);

  // if(debug) ShowMatrix(Matrix,n,4);
  if(debug) printf("Hole Solver...\n");

  if(type==0) tridag(Matrix[1],Matrix[2],Matrix[3],Matrix[4],v,n);
  else if(type==1) {
    double ** a = matrix(1,n,1,3);
    double ** al = matrix(1,n,1,1);
    int * indx = ivector(1,n);
    double d;
    for(i=1;i<=n;i++) for(j=1;j<=3;j++) a[i][j] = Matrix[i][j];
    for(i=0;i<n;i++) v[i] = Matrix[i][4];
    bandec(a,n,1,1,al,indx,&d);
    banbks(a,n,1,1,al,indx,v);
    free_ivector(indx,1,n);
    free_matrix(a,1,n,1,3);
    free_matrix(al,1,n,1,1);
  }
  else if (type==2) {
    i=0;
    for(i=1;i<=n;i++) v[i] = exp(phip(i+1)/VT);
    v[1] = (Matrix[1][4]-Matrix[1][3]*v[2])/Matrix[1][2];
    for(i=2;i<n;i++) v[i]= (Matrix[i][4] -Matrix[i][1]*v[i-1]-Matrix[i][3]*v[i+1])/Matrix[i][2];
    v[n] = (Matrix[n][4]-Matrix[n][1]*v[n-1])/Matrix[n][2];
  }
  else nrerror("Select right hole solver type!");
  
  // if(debug) ShowVector(v,n);

  error = 0.0;
  for(i=1;i<=n;i++) {
    tmp = fabs(phip(i+1)-VT*log(v[i]));
    if(tmp>1e3) {
      printf("%d: phip=%le  new=%le\n",i,phip(i+1),VT*log(v[i]));
    }
    error += tmp;
    set_phip(i+1,(1.0-x)*phip(i+1) + x*VT*log(v[i]));
  }

  free_vector(v,1,n);
  free_matrix(Matrix,1,n,1,4);

  if(debug) printf("Hole Solver done\n");

  return error;

}

/* -------------------------------------------------- */

/* Convert notation (k-layer,i-grid) into notation 1-xnum_mesh_points  */
 int Device::LayerToMesh(int kk, int ii) {
  int k,tmp =0;
  for(k=0;k<kk;k++) tmp += lyr[k].numx();
  tmp += ii+1;
  return tmp;
}

/* -------------------------------------------------- */

void Device::MeshToLayer(int j, int & k, int & i) {
  int tmp = j;
  k=-1;
  do {
    k++;
    tmp -= lyr[k].numx();
  } while(tmp>0);
  tmp += lyr[k].numx();
  i = tmp-1;
}

/* -------------------------------------------------- */

double Device::n(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].n();
}

/* -------------------------------------------------- */

double Device::p(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].p();
}

/* -------------------------------------------------- */

double Device::ni(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].ni();
}

/* -------------------------------------------------- */

double Device::n1(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].n1();
}

/* -------------------------------------------------- */
 
double Device::p1(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].p1();
}
/* -------------------------------------------------- */
 
double Device::phi(int num) {
   int k,i;
   MeshToLayer(num,k,i);
   return lyr[k].grid1d[i].ret_phi();
}
 
/* It goes from 1-to-xnum_mesh_points or 1-to-ynum_mesh_points */
double Device::phi2D(int numx, int numy) {
   int k,i;
   MeshToLayer(numx,k,i);
   if(numy<1 || numy>ynum_mesh_points) nrerror("PHI(): Y-coordinate out of range!");
   return lyr[k].grid2d[i][numy-1].ret_phi();
 }

void Device::set_phi(int numx, int j, double val) {
   int k,i;
  MeshToLayer(numx,k,i);
  lyr[k].grid2d[i][j-1].set_phi(val);
}

void Device::set_phi(int num, double val) {
  int k,i;
  MeshToLayer(num,k,i);
  lyr[k].grid1d[i].set_phi(val);
}

/* -------------------------------------------------- */

double Device::phin(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].ret_phi_n();
}

void Device::set_phin(int num, double val) {
  int k,i;
  MeshToLayer(num,k,i);
  lyr[k].grid1d[i].set_phi_n(val);
}

double Device::phin(int numx, int j) {
  int k,i;
  MeshToLayer(numx,k,i);
  return lyr[k].grid2d[i][j-1].ret_phi_n();
}

void Device::set_phin(int numx ,int j, double val) {
  int k,i;
  MeshToLayer(numx,k,i);
  lyr[k].grid2d[i][j-1].set_phi_n(val);
}

/* -------------------------------------------------- */

double Device::nf(int numx, int numy, double xf) {

  int k,i,j;

   MeshToLayer(numx,k,i);
   j = numy-1;
   return lyr[k].grid2d[i][j].nf(xf);

}

double Device::dnf(int numx, int numy, double xf) {

  int k,i,j;

   MeshToLayer(numx,k,i);
   j = numy-1;
   return lyr[k].grid2d[i][j].dnf(xf);

}

double Device::pf(int numx, int numy, double xf) {

  int k,i,j;

   MeshToLayer(numx,k,i);
   j = numy-1;
   return lyr[k].grid2d[i][j].pf(xf);

}

double Device::dpf(int numx, int numy, double xf) {

  int k,i,j;

   MeshToLayer(numx,k,i);
   j = numy-1;
   return lyr[k].grid2d[i][j].dpf(xf);

}


/* -------------------------------------------------- */

double Device::phip(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].ret_phi_p();
}

void Device::set_phip(int num, double val) {
  int k,i;
  MeshToLayer(num,k,i);
  lyr[k].grid1d[i].set_phi_p(val);
}

double Device::phip(int numx, int j) {
  int k,i;
  MeshToLayer(numx,k,i);
  return lyr[k].grid2d[i][j-1].ret_phi_p();
}

void Device::set_phip(int num, int j, double val) {
  int k,i;
  MeshToLayer(num,k,i);
  lyr[k].grid2d[i][j-1].set_phi_p(val);
}

/* -------------------------------------------------- */

double Device::Jn(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].ret_Jn();
}

void Device::set_Jn(int num, double val) {
  int k,i;
  MeshToLayer(num,k,i);
  lyr[k].grid1d[i].set_Jn(val);
}

/* -------------------------------------------------- */

double Device::Jp(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].ret_Jp();
}

void Device::set_Jp(int num, double val) {
  int k,i;
  MeshToLayer(num,k,i);
  lyr[k].grid1d[i].set_Jp(val);
}

/* -------------------------------------------------- */

double Device::D(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].ret_D();
}

void Device::set_D(int num, double val) {
  int k,i;
  MeshToLayer(num,k,i);
  lyr[k].grid1d[i].set_D(val);
}
/* -------------------------------------------------- */

double Device::SHR_Rec(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].SRH_Rec();
}

double Device::Trap_Rec_e(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].Trap_Rec_e();
}

double Device::Trap_Rec_h(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].Trap_Rec_h();
}

double Device::Opt_Rec(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].Opt_Rec();
}

double Device::QW_Rec(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].QW_Rec();
}

double Device::doping(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].ret_doping();
}

double Device::eps(int numx, int numy) {
  int k,i;
  MeshToLayer(numx,k,i);
  return lyr[k].grid2d[i][numy-1].eps();
}

double Device::chi(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].chi();
}

double Device::chi(int numx, int numy) {
  int k,i,j;
  MeshToLayer(numx,k,i);
  j = numy-1;
  return lyr[k].grid2d[i][j].chi();
}

double Device::Eg(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].Eg();
}

double Device::Nc(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].Nc();
}

double Device::Nv(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].Nv();
}

double Device::mobn(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].mob_n();
}

double Device::mobp(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].mob_p();
}
double Device::Dn(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].Dn();
}

double Device::Dp(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].Dp();
}

double Device::Tn(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].Tn();
}

double Device::Tp(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].Tp();
}

double Device::Tntrap(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].agrid->Tntrap();
}

double Device::Tptrap(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].agrid->Tptrap();
}

double Device::me(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].me();
}

double Device::mh(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].mh();
}


double Device::dx(int num) {
  int k,i;
  MeshToLayer(num,k,i);
  return lyr[k].grid1d[i].xsize();
}

double Device::dx(int numx, int numy) {
  int k,i;
  MeshToLayer(numx,k,i);
  return lyr[k].grid2d[i][numy-1].xsize();
}

double Device::dy(int numx, int numy) {
  int k,i;
  MeshToLayer(numx,k,i);
  return lyr[k].grid2d[i][numy-1].ysize();
}

double Device::Polar(int numx, int numy) {
  int k,i;
  MeshToLayer(numx,k,i);
  return lyr[k].grid2d[i][numy-1].Polarization();
}

double Device::mobn(int numx, int numy) {
  int k,i,j;
  MeshToLayer(numx,k,i);
  j = numy-1;
  return lyr[k].grid2d[i][j].mob_n();
}

double Device::mobp(int numx, int numy) {
  int k,i,j;
  MeshToLayer(numx,k,i);
  j = numy-1;
  return lyr[k].grid2d[i][j].mob_p();
}

double Device::dy(int num) {
  return lyr[0].ysize()/(double)ynum_mesh_points;
}

double Device::dt(int num) {
  return lyr[0].tsize()/(double)tnum_mesh_points;
}

/* --------------------------------------------------------------- */

double ** Device::ElectronDiscretization1D_2(short type, int & num) {

  int i;
  double a1,a2,a3,coeff;
  double ** Matrix;
  
  /* ALLOCATE MATRIX */

  num = xnum_mesh_points-2;
  Matrix = matrix(1,num,1,4);

  for(i=2;i<=num+1;i++) {

    /* left current coefficient */
    a1 = (phi(i)-phi(i-1))/VT + (chi(i)-chi(i-1))/VT + log(Nc(i)/Nc(i-1));
    a2  = fabs(a1) < 1e-10 ? 1 : a1 / (exp(a1) - 1);
    a3	= 2.0*a2*Dn(i)*Nc(i)*exp((phi(i)+chi(i))/VT);
    Matrix[i-1][1]= 2.0*a3/(dx(i)+dx(i-1));

    /* right current coefficient */
    a1 = ( phi(i)- phi(i+1))/VT + (chi(i)-chi(i+1))/VT + log(Nc(i)/Nc(i+1));
    a2  = fabs(a1) < 1e-10 ? 1 : a1 / (exp(a1) - 1);
    a3	= 2.0*a2*Dn(i)*Nc(i)*exp((phi(i)+chi(i))/VT);
    Matrix[i-1][3]= 2.0*a3/(dx(i)+dx(i+1));

    /* middle coefficient */
    Matrix[i-1][2] = - Matrix[i-1][1] - Matrix[i-1][3];
    coeff = (  ni(i)*ni(i) / ( Tn(i)*(p(i)+ni(i))+ Tp(i)*(n(i)+ni(i))) )*(2.0*dx(i)+dx(i-1)+dx(i+1));
    Matrix[i-1][2] -= exp(phip(i)/VT)*coeff;

    /* Known coefficient */
    Matrix[i-1][4] =  coeff;
  }

  /* LEFT BOUNDARY */

  Matrix[1][4] -= Matrix[1][1]*exp(-phin(1)/VT);
  Matrix[1][1] = 0;

  /* RIGHT BOUNDARY */

  Matrix[num][4] -= Matrix[num][3]*exp(-phin(num+2)/VT);
  Matrix[num][3] = 0;

  return Matrix;
}


/* --------------------------------------------------------------- */

double ** Device::HoleDiscretization1D_2(short type, int & num) {

  int i;
  double a1,a2,a3,coeff;
  double ** Matrix;
  

  /* ALLOCATE MATRIX */

  num = xnum_mesh_points-2;
  Matrix = matrix(1,num,1,4);

  for(i=2;i<=num+1;i++) {

    /* left current coefficient */
    a1 = (-phi(i)+phi(i-1))/VT + (-chi(i)+chi(i-1))/VT + (-Eg(i)+Eg(i-1))/VT + log(Nv(i)/Nv(i-1));
    a2 = fabs(a1) < 1e-10 ? 1.0 : a1 / (exp(a1) - 1.0);
    a3 = 2.0*a2*Dp(i)*Nv(i)*exp(-(phi(i)+chi(i)+Eg(i))/VT);
    Matrix[i-1][1]= 2.0*a3/(dx(i)+dx(i-1));

    /* right current coefficient */
    a1 = (-phi(i)+phi(i+1))/VT + (-chi(i)+chi(i+1))/VT + (-Eg(i)+Eg(i+1))/VT + log(Nv(i)/Nv(i+1));
    a2  = fabs(a1) < 1e-10 ? 1.0 : a1 / (exp(a1) - 1.0);
    a3 = 2.0*a2*Dp(i)*Nv(i)*exp(-(phi(i)+chi(i)+Eg(i))/VT);
    Matrix[i-1][3]= 2.0*a3/(dx(i)+dx(i+1));

    /* middle coefficient */
    Matrix[i-1][2] = - Matrix[i-1][1] - Matrix[i-1][3];
    coeff = ( ni(i)*ni(i) / ( Tn(i)*(p(i)+ni(i))+ Tp(i)*(n(i)+ni(i))) )*(2.0*dx(i)+dx(i-1)+dx(i+1));
    Matrix[i-1][2] -= exp(-phin(i)/VT)*coeff;  

    /* known coefficient */
    Matrix[i-1][4] = coeff;

  }

  /* LEFT BOUNDARY */

  Matrix[1][4] -= Matrix[1][1]*exp(phip(1)/VT);
  Matrix[1][1] = 0;    

  /* RIGHT BOUNDARY */

  Matrix[num][4] -= Matrix[num][3]*exp(phip(num+2)/VT);
  Matrix[num][3] = 0;

  return Matrix;
}


/* --------------------------------------------------------------- */
/* --------------------------------------------------------------- */

void Device::Current(short type, short debug) {

  int i;

  if(debug) printf("Current Solver...\n");

  if(type==0) {
    if(debug) printf("METHOD 1: the quasi-fermi levels.\n");
    if(debug) printf("This method supposes that if one boundary is n-doped the other is p-doped\n");
    
    if(doping(1)>0.0) {
      set_Jn(1,0.0);
      for(i=2;i<=xnum_mesh_points;i++) set_Jn(i,-2.0*Q*n(i)*mobn(i)*(phin(i)-phin(i-1))/(dx(i)+dx(i-1)));
      set_Jp(xnum_mesh_points,0.0);
      for(i=xnum_mesh_points-1;i>=1;i--) set_Jp(i,-2.0*Q*p(i)*mobp(i)*(phip(i+1)-phip(i))/(dx(i)+dx(i+1)));
    }
    else {
      set_Jn(xnum_mesh_points,0.0);
      for(i=xnum_mesh_points-1;i>=1;i--) set_Jn(i,-2.0*Q*n(i)*mobn(i)*(phin(i+1)-phin(i))/(dx(i)+dx(i+1)));
      set_Jp(1,0.0);
      for(i=2;i<=xnum_mesh_points;i++) set_Jp(i,-2.0*Q*p(i)*mobp(i)*(phip(i)-phip(i-1))/(dx(i)+dx(i-1)));
    }
  }
  else if(type==1) {
    if(debug) printf("METHOD 2: recombination SRH\n");
    if(sim->traps) {
      if(doping(1)>0.0) {
	set_Jn(1,0.0);
	for(i=2;i<=xnum_mesh_points;i++) set_Jn(i,Jn(i-1)+Q*0.5*(dx(i)+dx(i-1))*(SHR_Rec(i)+Trap_Rec_e(i)+Trap_Rec_h(i)+Opt_Rec(i)));
	set_Jp(xnum_mesh_points,0.0);
	for(i=xnum_mesh_points-1;i>=1;i--) set_Jp(i,Jp(i+1)+Q*0.5*(dx(i)+dx(i+1))*(SHR_Rec(i)+Trap_Rec_e(i)+Trap_Rec_h(i)+Opt_Rec(i)));
      }
      else {
	set_Jp(1,0.0);
	for(i=2;i<=xnum_mesh_points;i++) set_Jp(i,Jp(i-1)-Q*0.5*(dx(i)+dx(i-1))*(SHR_Rec(i)+Trap_Rec_e(i)+Trap_Rec_h(i)+Opt_Rec(i)));
	set_Jn(xnum_mesh_points,0.0);
	for(i=xnum_mesh_points-1;i>=1;i--) set_Jn(i,Jn(i+1)-Q*0.5*(dx(i)+dx(i+1))*(SHR_Rec(i)+Trap_Rec_e(i)+Trap_Rec_h(i)+Opt_Rec(i)));
      }
    }
    else {
      if(doping(1)>0.0) {
	set_Jn(1,0.0);
	for(i=2;i<=xnum_mesh_points;i++) {
	  set_Jn(i,Jn(i-1)+Q*0.5*(dx(i)+dx(i-1))*(SHR_Rec(i)+Opt_Rec(i)));
	}
	set_Jp(xnum_mesh_points,0.0);
	for(i=xnum_mesh_points-1;i>=1;i--) {
	  set_Jp(i,Jp(i+1)+Q*0.5*(dx(i)+dx(i+1))*(SHR_Rec(i)+Opt_Rec(i)));
	}
      }
      else {
	set_Jp(1,0.0);
	for(i=2;i<=xnum_mesh_points;i++) set_Jp(i,Jp(i-1)-Q*0.5*(dx(i)+dx(i-1))*(SHR_Rec(i)+Opt_Rec(i)));
	set_Jn(xnum_mesh_points,0.0);
	for(i=xnum_mesh_points-1;i>=1;i--) set_Jn(i,Jn(i+1)-Q*0.5*(dx(i)+dx(i+1))*(SHR_Rec(i)+Opt_Rec(i)));
      }
    }
  }
  else nrerror("Select a method for current");

  if(debug) printf("Current Solver done\n");

}

/* --------------------------------------------------------------- */
/* --------------------------------------------------------------- */

double Device::GrabCurrent(short debug) {

  int i;
  double tmp = 0.0;

  for(i=1;i<=xnum_mesh_points;i++) tmp += Jn(i)+Jp(i);
  tmp = tmp/xnum_mesh_points;
  return tmp;

}

/* --------------------------------------------------------------- */
/* --------------------------------------------------------------- */

int Device::simulator() {

  int count,i;
  char  ch=' ';
  double e_err,h_err,bias_tmp, * bias_dstr;
  FILE * fp = 0;

  /* DUMP INITIAL CONDITION FILES */
  dump_files('0');

  if(sim->traps) printf("WITH TRAPS\n");

  /* CHECK IF THE BIAS SETTING IS CONSISTENT */
  if(sim->final_bias!=sim->initial_bias) {
    printf("Initial bias = %lf, final bias = %lf, bias step = %lf\n",sim->initial_bias,sim->final_bias,sim->step_bias);
    if((sim->final_bias>sim->initial_bias && sim->step_bias<0) || (sim->final_bias<sim->initial_bias && sim->step_bias>0) ) {
      printf("Final and initial bias setting contradict the sign of the bias step\n Changing sign ...\n");
      printf("New setting: initial bias = %lf, final bias = %lf, bias step = %lf\n",sim->initial_bias,sim->final_bias,-sim->step_bias);
      sim->step_bias = -sim->step_bias;
      printf("Are you happy with that?[y/n] ");
      while(ch!='y' && ch!='n') scanf_s("%c",&ch);
      if(ch=='n') nrerror("Abort");
    }
    fopen_s(&fp,"OUTPUT_FILES/iv_curve.dat","w");
    fclose(fp);
  }
  else sim->step_bias = 1;
  
  bias_tmp = sim->initial_bias;


  while(bias_tmp<=sim->final_bias) {

    printf("\n\nBIAS = %lf\n\n",bias_tmp);
    count = 0;
    do {
      count++;
      if(ynum_mesh_points<=1) { // 1-D Simulation
		if(count%sim->num_poisson==0) 
			Poisson_solvde(sim->debug);
		if(count%sim->num_schrodinger==0) 
			Schrodinger(sim->debug);
		if(count%sim->num_continuity==0)  {
			e_err = ElectronSolver(2,1.0,sim->debug);
			h_err = HoleSolver(2,1.0,sim->debug);
		}

		if(count%sim->num_current==0) Current(1,sim->debug);
		if(count%sim->num_error==0) printf("%d\tn=%le\t p=%le\n",count,e_err,h_err);
		if(count%sim->num_file==0) dump_files('x');
      }
      else if(ynum_mesh_points>1) { // 2-D Simulation
		if(count%sim->num_poisson==0) 
			Poisson_2D();
		if(count%sim->num_continuity==0) 
			Continuity_2D();
      }
      else nrerror("DEVICE:simulator No dimension specified!");

    } while(count<sim->num_iteration && (e_err>sim->precision || h_err>sim->precision));


    fopen_s(&fp,"OUTPUT_FILES/iv_curve.dat","a");
    fprintf(fp,"%lf\t%le\n",bias_tmp,GrabCurrent(sim->debug));
    fclose(fp);
    bias_tmp += sim->step_bias;
    if(bias_tmp<=sim->final_bias) {
      bias_dstr = bias_distribution(bias_tmp,sim->debug);
      for(i=0;i<num_lyr;i++) lyr[i].bias_modify(bias_dstr[i+1]);
      free_vector(bias_dstr,1,num_lyr);
    }

  }
  
  dump_files('f');

  return 0;  

}

// --------------------------------------------------------------- 
// ---------------------------------------------------------------

// This routine sets the difference of potential among the
// existing layers according to the doping of the layer and the 
// applied bias. It consider the fact that if there is not a 
// large doping difference, the potential has to be the same.
// The bias is supposed to be applied to the end of the device
// that is the last layer in the file device.inp. The bias is
// applied half on the right side and the other half in the right one.

double * Device::bias_distribution(double bias_inp, short debug) {

  int i;
  double * bias = vector(1,num_lyr);

  bias[1] = - 0.5 * bias_inp;
  bias[num_lyr] =  0.5 * bias_inp;

  for(i=1;i<num_lyr-1;i++) {
    if(lyr[0].ret_doping_rhs()>0 && lyr[i].ret_doping_lhs()>0) {
      if(lyr[i].ret_doping_lhs()*1e2<lyr[0].ret_doping_rhs()) bias[i+1] = 0;
      else bias[i+1] = bias [1];
    }
    else if(lyr[0].ret_doping_rhs()<0 && lyr[i].ret_doping_lhs()<0) {
      if(-lyr[i].ret_doping_lhs()*1e2<-lyr[0].ret_doping_rhs()) bias[i+1] = 0;
      else bias[i+1] = bias [1];
    }
    else if(lyr[num_lyr-1].ret_doping_lhs()<0 && lyr[i].ret_doping_rhs()<0) {
      if(-lyr[i].ret_doping_rhs()*1e2<-lyr[num_lyr-1].ret_doping_lhs()) bias[i+1] = 0;
      else bias[i+1] = bias [num_lyr];
    }
    else if(lyr[num_lyr-1].ret_doping_lhs()>0 && lyr[i].ret_doping_rhs()>0) {
      if(lyr[i].ret_doping_rhs()*1e2<lyr[num_lyr-1].ret_doping_lhs()) bias[i+1] = 0;
      else bias[i+1] = bias [num_lyr];
    }
    else {
      printf("%d = %le, %le, %le\n",i,lyr[0].ret_doping_rhs(),lyr[i].ret_doping_lhs(),lyr[num_lyr-1].ret_doping_lhs());
      nrerror("Something wrong with the bias setting!");
    }
  }
  
  return bias;

}


/* ----------------------------------------------------- */
/* ----------------------------------------------------- */


